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ABSTRACT 

Aims. We apply various expansion schemes that may be used to study gravitational clustering to the simple case ol the 
Zeldovich dynamics. 

Methods. Using the well-known exact solution of the Zeldovich dynamics we can compare the predictions of these 
various perturbative methods with the exact nonlinear result. We can also study their convergence properties and their 
behavior at high orders. 

Results. We find that most systematic expansions fail to recover the decay of the response function in the highly nonlinear 
regime. "Linear methods" lead to increasingly fast growth in the nonlinear regime for higher orders, except for Pade 
approximants that give a bounded response at any order. "NonUnear methods" manage to obtain some damping at 
one-loop order but they fail at higher orders. Although it recovers the exact Gaussian damping, a resummation in the 
high-fc limit is not justified very well as the generation of nonlinear power does not originate from a finite range of 
wavenumbers (hence there is no simple separation of scales). No method is able to recover the relaxation of the matter 
power spectrum on highly nonlinear scales. It is possible to impose a Gaussian cutoff in a somewhat ad-hoc fashion 
to reproduce the behavior of the exact two-point functions for two different times. However, this cutoff is not directly 
related to the clustering of matter and disappears in exact equal-time statistics such as the matter power spectrum. 
On a quantitative level, on weakly nonlinear scales, the usual perturbation theory, and the nonlinear scheme to which 
one adds an ansatz for the response function with such a Gaussian cutoff, are the two most efficient methods. We can 
expect these results to hold for the gravitational dynamics as well (this has been explicitly checked at one- loop order), 
since the structure of the equations of motion is identical for both dynamics. 
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1. Introduction 

■ The growth of large-scale structures in the Universe 
through the amplification of small primordial fluctuations 

. by gravitational instability is a key ingredient of mod- 

■ ern cosmology (Peebles 1980). This process can be used 
, to constrain cosmological parameters through the depen- 
dence of the matter power spectrum on scale and redshift. 
For instance, observations of highly nonlinear objects such 

, as galaxy clusters can help constrain the normalization of 
the matter power spectrum and the average matter density 
(Oukbir & Blanchard 1992; Younger et al. 2005). Although 
theoretical predictions are not very accurate on these small 
scales, one can still derive useful constraints because they 
are very rare objects, so that their dependence on cosmo- 
logical parameters is very steep. Alternative probes, such as 
baryonic acoustic oscillations (Eisenstein et al. 1998, 2005) 
or weak lensing surveys (Munshi et al. 2007; Massey et al. 
2007), focus on weakly nonlinear scales where they aim at 
measuring the matter distribution through its power spec- 
trum or its two-point correlation, which are sensitive to typ- 
ical fluctuations. In such cases, one needs an accurate theo- 
retical prediction to derive useful constraints on cosmology. 
This problem is usually tackled through N-body simula- 
tions or perturbative expansions for scales not too far from 
the linear regime (Bernardeau et al. 2002). However, nu- 
merical simulations are rather costly and analytical meth- 
ods may have the advantage of leading to a better under- 
standing of the physics at work. Moreover, on such large 



scales, a hydrodynamical description should be sufficient 
(i.e. one neglects shell crossing). This facilitates analytical 
approaches as one can use the continuity and Euler equa- 
tions instead of the Vlasov equation. Therefore, there has 
recently been renewed interest in perturbation theory to 
devise analytical methods that could exhibit good accu- 
racy on weakly nonlinear scales to be used for such obser- 
vational probes (Crocce & Scoccimarro 2006a, b; Valageas 
2007; McDonald 2007; Matarrese & Pietroni 2007a,b). 

Thus, Crocce & Scoccimarro (2006a, b) find that one 
can perform a partial resummation of the diagrammatic 
series that appears in the standard perturbative expansion 
to obtain a response function that decays into the nonlinear 
regime as expected (whereas the standard expansion grows 
as a polynomial of increasing order as we truncate the se- 
ries at higher order). Moreover, this result agrees well with 
numerical simulations and can be used as an intermediate 
tool for obtaining a more accurate prediction of the matter 
power spectrum than with the usual perturbation theory 
(Crocce & Scoccimarro 2007). On the other hand, Valageas 
(2007) present a path-integral formalism, so that the system 
is fully defined by its action S (or its partition function Z). 
Then, one can apply the usual tools of field theory, such 
as large- A'^ expansions (similar to a semi-classical expan- 
sion over powers of 7i or a generalization to N fields), to 
derive the matter power spectrum (Valageas 2007). Note 
that this method also applies to the highly nonlinear scales 
described by the Vlasov equation (Valageas 2004). Next, 
Matarrese & Pietroni (2007a) have recently proposed an 
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alternative method based on the path-integral formalism 
where one considers the dependence of the system on a 
large-wavenumber cutoff A. This gives rise to a new set of 
equations and, by taking the limit A ^ oo, one recovers the 
original system. These various methods may be seen as dif- 
ferent reorganizations of the standard perturbative expan- 
sion. Of course, they all involve some truncation at some 
stage (otherwise the problem would be solved exactly), and 
they are all consistent up to that order (i.e. they only differ 
by higher-order terms). 

To check the range of validity of such expansions, one 
must compare their predictions with N-body simulations 
(and assume for observational purposes that the accuracy 
remains the same for close cosmological parameters). This 
is not very convenient, since simulations themselves may be 
of limited accuracy. Along with this problem, the behav- 
iors of other two-point functions than the power spectrum 
P(k), such as the response function and self-energies and es- 
pecially different-time functions such as ((5(xi, ti)(5(x2, i2)) 
with ti 7^ ^2, have not been analyzed in detail from N- 
body simulations (which do not give direct access to self- 
energies either). Therefore, it is interesting to investigate 
these theoretical methods applied to a simpler dynamics 
that can be solved exactly. Then, one can do a detailed 
comparison of the predictions of various expansion schemes 
with the exact nonlinear results. Moreover, one can recon- 
struct such expansion schemes in a direct manner from the 
exact two-point functions without computing high-order 
diagrams that involve many integrals, simply by expand- 
ing back the exact nonlinear result. In this way one can 
more easily investigate the convergence properties of these 
expansions and the behavior of high-order terms. A sim- 
ple dynamics that can be solved exactly but that remains 
close to the gravitational dynamics (at least up to weakly 
nonlinear scales) is provided by the Zeldovich approxima- 
tion (Zeldovich 1970; Gurbatov et al. 1989; Shandarin & 
Zeldovich 1989). The latter was originally devised as an 
approximation to the gravitational dynamics. Here we take 
a different point of view as we modify the equations of 
motion so that the system is exactly given by the simple 
Zeldovich dynamics. Then, we apply to these new equa- 
tions of motion various methods which can be applied to 
both dynamics (and to other stochastic dynamics such as 
the Navier-Stokes equations). Taking advantage of the ex- 
act results which can be obtained for the Zeldovich dynam- 
ics and its simpler properties we study the accuracy and 
the general properties of these expansion schemes in detail. 
This should shed light on the behavior of these methods ap- 
plied to the gravitational dynamics, because both dynamics 
exhibit similar equations of motion, and these expansions 
apply in identical manner to both systems. 

This article is organized as follows. First, in Sect. [2] 
we derive the equations of motion associated with the 
Zeldovich dynamics and their linear solution. Next, in 
Sect. [3] we obtain the path-integral formulation of this sys- 
tem, starting from either the differential form or the in- 
tegral form of the equations of motion, in order to make 
the connection with the different approaches used in the 
literature. Then, we briefly describe how some expansion 
schemes can be built from this path-integral formalism, 
such as the large- iV expansions in Sect. [3] and the evolu- 
tion equations with respect to a high-fc cutoff in Sect. [S] 
Before investigating such methods, we first derive the exact 
nonlinear two-point functions which can be obtained from 



the well-known exact solution of the Zeldovich dynamics 
in Sect. [S] Then, we describe the behavior of the standard 
perturbation theory in Sect. [7] and of the steepest-descent 
method (built from a large- TV expansion) in Sect. [S] Next, 
we discuss in Sect. 1^] the high-/c resummation proposed by 
Crocce & Scoccimarro (2006b) to improve the behavior of 
such expansion schemes. We turn to the 2PI effective action 
method in Sect. [TO] (a second approach built from a large- iV 
expansion) and to simple nonlinear schemes associated with 
this expansion in Sect. [TTl We investigate simple nonlinear 
schemes associated with the dependence on a high-fc cutoff 
in Sect. [12] Finally, in Sect. [13] we study the quantitative 
predictions on weakly nonlinear scales of these methods at 
one- loop order and we conclude in Sect. 1141 

2. Equations of motion 

2.1. Zeldovich approximation 

On scales much larger than the Jeans length, both the cold 
dark matter and the baryons can be described as a pres- 
sureless dust. Then, we can neglect orbit crossings and use 
a hydrodynamical description governed by the equations of 
motion (Peebles 1980): 



|^+V.[(H-<5)v]=0, 



9v 

a7 



7i:v+ (v.V)v 



-V(/., 



A0 = -n^n's, 



(1) 



(2) 



(3) 



where t — Jdt/a is the conformal time (and a the scale 
factor), Ti = dlna/dr the conformal expansion rate, and 
fim the matter density cosmological parameter. Here, 5 is 
the matter density contrast and v the peculiar velocity. 
Since the vorticity field decays within linear theory (Peebles 
1980), we take the velocity to be a potential field so that 
V is fully specified by its divergence 9 or by its potential x 
with 



— V.v, V = — Vx whence 6 = —Ax. 



(4) 



In the linear regime, one finds that the linear growing mode 
satisfies 



h = -fH^L whence (j)L = ^^"f^ XL, 



(5) 



where /(r) is defined from the linear growing rate Dj^{t) 
of the density contrast by 



d\nD+ 1 dlnDi 



/ = _ 

din a 7Y dr ' 

and (r) is the growing solution of 



dr2 



n- 



.dD: 



dr 



(6) 



(7) 



If we make the approximation that relation ([5]) remains 
valid in the nonlinear regime, that is, we replace the Poisson 
equation Q by the second Eq.®: (f> = SflynTi-x/^f , then 
we obtain for the Euler equation ([2]): 



(8) 
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Obviously, as shown by Eq.®, within this approximation 
the velocity field now evolves independently of the density 
field. As is well known (Gurbatov et al. 1989), approxi- 
mation ([8| is actually identical to the Zeldovich approxi- 
mation. Indeed, a change of variables for the velocity field 
yields 



du 



+ (u.V)u = with V 



fdD^ 



\ dr 



(9) 



Equation ^ is the equation of motion of free particles, 
du/dZ?-|_ ~ 0, hence the trajectories are given by 



X = q + -D+(T)sLo(q), V 



~d7 



SLo(q), 



(10) 



where q is the Lagrangian coordinate and 8 = 8^ = D^Si^q 
is the displacement field that is exactly given by the linear 
theory. Equation (fTO|) is the usual definition of the Zeldovich 
approximation (i.e. setting s = s^). 

Thus, the Zeldovich approximation corresponds to a 
change in the linear term of the Euler equation, keeping 
the quadratic term and the continuity equation unchanged. 
Therefore, the analysis presented in Valageas (2007) for 
the case of the exact gravitational dynamics applies to the 
Zeldovich dynamics up to minor modifications. First, the 
equations of motion ^ and (O read in Fourier space as 



In Eq. (fT7)) and in the following, we use the convention that 
repeated coordinates are integrated over as 

0{x, x').i;{x') = / dk'dry' ^ O,,,' (k, rj; k', r;')^»' (k', 7/). (18) 

i' — l 

The matrix O reads 

0(x,x')=(^f o^\^6D{k-k')SD{Ti-v') (19) 

whereas the symmetric vertex Ks(x; xi,X2) — Ks(x; X2,xi) 
writes as 

Ks{x\xi,X2) = ^^(ki -f k2 - k)5D(?7i - ?y)<^D(?72 - 77) 

X7^^,,^(ki,k2) (20) 

with 

, , X «(k2,ki) «(ki,k2) 
7i;i,2(iii7k2) = Ty > 7i;2,i(ki,k2j = , (21) 



72;2.2(kl,k2)=/3(ki,k2), 



(22) 



9(5(k, r) 
9t 



^(k, ^) = - J dkidk2 Sd(Mi + k2 - k) 

xa(ki,k2)0(ki,T)(5(k2,r) (11) 



-J dkidk2 ^^(ki + k2 - k)/3(ki, k2)0(ki, r)0(k2, r) (12) 

where is the Dirac distribution. The coupling functions 
a and /3 are given by 

(ki+Mki ^ki+k2p(ki^ 
a(ki,k2) = 72 ,P(ki,k2) = —j-^ ,(13) 



and zero otherwise (Crocce & Scoccimarro 2006a). We can 
note that all the dependence on cosmology is contained in 
the time-redshift relation 77 z. Indeed, the equation of 
motion (|17p written in terms of the coordinate 77 no longer 
involves time-dependent factors such as flm/ f^- Therefore, 
the evolution of the density field only depends on cosmol- 
ogy through the time-coordinate r]{z). In this article we 
study the system defined by the the equation of motion 
((T7]), which shows the exact solution PH)) . This will allow 
us to compare various expansion methods with exact non- 
linear results. 



2.2. Linear regime 

On large scales or at early times where the density and ve- 
locity fluctuations are small, one can linearize the equation 
of motion P7|) . which yields O.tpL = 0. This gives the two 
linear modes: 



and we defined the Fourier transforms as 
dx 



(14) 



^+ = e" ( : 1 , = 



(23) 



As in Crocce & Scoccimarro (2006a, b), let us define the 
two-component vector ip as 



V'(k, v) 



V'i(k, v) 
7^2 (k, 77) 



<5(k,77) 

-eik, v)/fn J ' 



(15) 



where we have introduced the time coordinate rj defined 
from the linear growing rate of the density contrast 
(normalized to unity today): 



77 = lni:>+(T) with D+{z = 0) = 1. (16) 
Then, the equations of motion (fTT|l - (fT^ can be written as = e''5io(k) 

0{x, x').i>{x') = Ks{x; Xl,X2)■^p{xl)^|;{x2), (17) 



Of course we recover the linear growing mode ip^ of the 
gravitational dynamics, since approximation ([5]) is valid in 
this case. However, the usual decaying mode ip- has been 
changed to a constant mode. As seen in Ea.(P5|). it corre- 
sponds to a mere perturbation of the density field that is 
transported by the unchanged velocity field. Indeed, since 
the velocity field is now decoupled from the density field, it 
obeys a first-order differential equation in the linear regime 
(rather than a second-order differential equation), which 
only admits one linear mode. As usual we define the initial 
conditions by the linear growing mode tpL- 



(24) 



where we have introduced the coordinate x = (k, 77, i) where 
I = 1, 2 is the discrete index of the two-component vectors. 



where SloO^) is the linear density contrast today at redshift 
z = 0. In this fashion the system pT|l - (ll2p that we study 
here agrees with the gravitational dynamics in the linear 
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regime. Besides, from Eqs.(IS]) and pi7|) . we see that the 
displacement field Sio(q) obeys 



Vq.SLo 



5L0- 



(25) 



Moreover, we assume Gaussian homogeneous and isotropic 
initial conditions defined by the linear power spectrum 
Pwik): 



('5L0(kl)^L0(k2)) = fe(ki + k2) PLo(fci)- 



(26) 



Then, as for the gravitational dynamics studied in Valageas 
(2007), the linear two-point correlation function Gl{xi, X2) 
reads as: 



Gl(xi,X2) = {'^L{xi)'il^L{x2)) 



<5z3(ki+k2)e"i+"^Pio(fci) ; { 



I 1 



(27) 



As in Valageas (2004), it is convenient to introduce the 
response function R{xi,X2) (related to the propagator used 
in Crocce & Scoccimarro (2006a, b), see Sect. 13.2"^ below) 
defined by the functional derivative 



R{Xi,X2) 



, H{xi) , 

'sax2) 



(28) 



where C(x) is a "noise" added to the r.h.s. in Eq. (fT7)l . Thus, 
R{xi, X2) measures the linear response of the system to an 
external source of noise. Because of causality it contains an 
Heaviside factor 9{rii — 772) since the field ip^xi) can only 
depend on the values of the "noise" at earlier times 772 < Vi ■ 
Moreover, it satisfies the initial condition: 



■q^ : i?(xi,X2) (5_D(ki - k2)(5ii,i 



(29) 



In the linear regime, the response Rl can be obtained from 
the initial condition (P^)) and the linear dynamics O.Rl — 
for 771 > 772 (as implied by the definition (pS)) and O.ipL — 
0). This yields (Crocce et al. 2006) 



Rl{xi,X2) = SoCki - k2)6'(77i - 772) 
l\ , / 1 -1 



1 







(30) 



This expression holds for any cosmology, whereas for the 
case of the gravitational dynamics factors, such as fim//^, 
lead to a small explicit dependence on cosmological param- 
eters. Note that by symmetry the two-point correlation G 
has the form 

G{xi,X2) Soi^i + k2)Gi^^i^{ki;r]i,r]2) (31) 
with 

Gii,t2{k;m.m) = ^'^^.^^(fc; 772, 771) (32) 

whereas the response function has the form 

R{xi,X2) = (5D(ki -k2)6'(77i - m) Rii.i2{h;m,V2)- (33) 

On the other hand, as noticed above, the linear two-point 
functions obey 

0{x,z).GL{z,y)^Q, 0{x,z).RL{z,y) = SD{x-y). (34) 

This can also be checked from the explicit expressions 
([?7|) . (|30p . Finally, it is convenient to define the power per 
logarithmic wavenumber A^(fc) by 

A2(fc) =47rfc3p(A:), A2(A:; 771, 772) = 4^fc3Gii(fc; 771, 772)(35) 



where the second expression generalizes A^(fc) at different 
times. Note that for 771 7^ 772 we can have A^(fc; 771, 772) < 0, 
whereas at equal times we have A^(k;r],ri) > 0. Then we 
have, for instance, 



k k\xi - X2I 



(36) 



('^(X1)<5(X2)) 

Thus, for a CDM cosmology the linear power Aj^^{k) grows 
as at low k and as In A; at high k. 

3. Path-integral formalism 

3.1. Differential form 

As in Valageas (2004, 2007) we can apply a path-integral 
approach to the hydrodynamical system described in 
Sect. 12.11 Let us briefiy recall how this can be done from 
the differential equation (fT7| (see also Martin et al. 1973; 
Phythian 1977). In order to explicitly include the initial 
conditions, we rewrite Eq.ljTT]) as 



with 7/^ = for 77 < 77/ and 
lii{x) = Soiv- 77/)e''^(5Lo(k) 



(37) 



SD{v-ViHi{x),m 



where we have introduced the coordinate x: 
x=(k,7) and ipjlx) ^ ipLix^rji). 



(39) 



Thus, the source m (which formally plays the role of some 
external noise) merely provides the initial conditions at 
time 77/, obtained from the linear growing mode (|24p . We 
shall eventually take the limit 77/ —00. Next, we define 
the generating functional Z[j] by 



(40) 



where we took the average over the Gaussian initial condi- 
tions 



{m) = 0, {lJ.l{xi)fIl{x2)) = Ai{xi,X2), 

with 



(41) 



A/(a:i,a:;2) = Soim - 'ni)SD{m - Vi)Giixi,X2), (42) 
Gi{xi,X2) = Gl(^i,77/;x2,77/). (43) 

All statistical properties of the field ip may be obtained 
from Z[j]. It is convenient to write Eq. (|10)) as 



Z[j] = J [d/ij] [dV-] I det M\ Sd(pi - O-V' + K,.i^i^) 

^gj.^-iM/.A7^Mi^ (44) 

where the Jacobian |detM| is defined by the functional 
derivative M = S^i/Sip. As in Valageas (2007), a simple 
computation shows that this Jacobian is equal to an irrele- 
vant constant. Then, introducing an imaginary ghost field 
A to express the Dirac as an exponential and performing 
the Gaussian integration over /i/, we obtain 



Z[j] = / [dV-lidA] e^-i^+^-(-o.4,+K^.4,4,)+^x.A,.x^ 



(45) 
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Thus, the statistical properties of the system pT)) are de- 
scribed by the action S[ijj, A] defined by 



Siijj, A] = X.{0.ip - Ks.H^) - -A.A/.A. 



(46) 



Moreover, we can note that adding a "noise" C, to the r.h.s. 
of Eq. (jl7p amounts to changing fij ^ fij + C, which trans- 
lates into 5 — !■ 5 — A-C- Therefore, functional derivatives 
with respect to C a-re equivalent to insertions of the ghost 
field A. In particular, we have 



X2) = (V'(a;i)A(x2)), (A) - 0, (AA) = 0. 



(47) 



The response function R is also related to the correlation 
with the initial conditions fij through 



S 



R.Ai 



Ar.X 



,\.i-0.rP+K,.TPrP) + ^X.Ai.X 



(48) 



since the integral of a total derivative vanishes, and we have 
used the symmetry of A/. This also reads as 



{ilj{xi)ipi{x2)) = R{xi;x,i]i) X Gi{x;x2) 



(49) 



where we define the cross-product x as the dot product 
IjlSp without integration over time, such as: 



Rx^i= /dk'^i?y(k,77;k' 



(50) 



3.2. Integral form 

In order to make the connection with the approach devel- 
oped in Crocce & Scoccimarro (2006a, b), we describe here 
how the same path-integral method can be applied to the 
equation of motion p7p written in integral form. 



3.2.1. Letting rji —oo 

First, as in Valageas (2001) (see also Scoccimarro 2000), we 
can integrate the equation of motion (jl7p as 



Ipix) = ll^hix) + Ks{x]Xi,X2).'tlj{xi)ljj{x2) 

with 

0.ks = Ks or K.^Rl.K, 



(52) 



as seen from Eq. ([M)) . Here the initial time rji no longer 
appears, because we have already taken the limit rjj ~oo. 
Then, following the same procedure as in Sect. 13.1] we can 
write 



Z[j] 



i.V[-0I.]-iVL.G-\^t 



(53) 



Introducing again an imaginary field x to impose the con- 
straint associated with the equation of motion (|51|1. we fi- 
nally obtain (the Jacobian is equal to unity): 



Tl) are now 



Thus, the statistical properties of the system 
described by the action iS['0,x] defined by 

S[i;,x]^X-{^~K,.iji;)-^X-GL-X- (55) 

Note that this formulation is equivalent to the one described 
in Sect. 13.11 except that we have already taken the limit 
r]i — oo directly into the equation of motion (|5ip . From 
the response field Xj we can again obtain a new response 
function TZ associated with Ea. ([5T|) . From the comparison 
with Eg. fiSl) . we see that we have the relations between 
both approaches: 

X-A.O, n^{yJx)^R-0^{^}, (56) 



where in the last expression we recall that from Eq. (|5ip a 
variation with respect to an external noise C can be seen 
as a variation with respect to tp^. In the linear regime we 
simply have TZL{x,y) — Sd{x ~ y). Moreover, in a fashion 
similar to Ea. (|49p we have the property 
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-^ + (^L-X 

ox 



X.{-rp+Ks.i.'tl') + ^X-GL.X 



= / m[dx]r^ 

J L "A J 

= mCL-x)) (57) 
which yields the relation 

{ip{xi)ipL{x2)) ^TZ{xi,x).Gl{x,X2), (58) 
where we use the symmetry oi G^. 

3.2.2. Integral form with finite r/j 

Finally, as in Crocce & Scoccimarro (2006a,b), it is possible 
to apply the initial conditions at some finite time 77/, as in 
Sect. 13.11 Thus, we may write the linear growing mode '0L 
at times 77 > 77/ as 



ipLix) = Rl{x;x','I]i) X tpi{x') 



(59) 

where ipj was defined in Eq. (|39p and the cross-product x in 
Eq. ([50| . Then, following the same procedure as in Eqs.fST])- 
([54)) . where the Gaussian average is now taken over the 
field ipi with two-point correlation G/, we now obtain the 
generating functional: 



(51) Z[?]=/[d7A/][dV'][dx]e^-'''+'^-(-^^'^'^^-'''+^=-''"^)-5'^^><GFV/ 



'dip] [dx] 



^3.4>+X-{~4'+Ks.4>4>) + ix-(RLXGixRl).x 



(60) 



Of course, we can check that, by taking the limit rjj —00 
in Eq. (pD|) . we recover Eq. (|Ml) since we have 

??i,'72>?7/: Gl(?7i, 772) = RL{m,Vi) x G/ x Rl{r]2,r]i).{61) 

The system is now described by the action S[ip,x] defined 
by 

X] = X-{^ - Ks.iPiP) - ^X-{Rl xGiX Rl).x. (62) 

Next, we can define a response function with respect to the 
initial conditions by 



(54) 'R{xi,X2) = { 



S^iix2)' 



{^{xi)x{x).Rl{x;x2,vi))- (63) 
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From the comparison of (p^ with we obtain x — 
and 

1i{xi,X2) = R.O.Rl = i?(xi;x2,77/)- (64) 

Thus, the response '^(a;i,X2), which is called the "prop- 
agator G'iji2(A:i,77i)(5£)(ki — k2)" in Crocce & Scoccimarro 
(2006a, b) is equal to the response function R of Sect. l3.1l re- 
strictcd to time 772 — Vii without taking the limit rjj — > ~oo. 
Finally, we can note that from Eq. (H5)) we have 

{ip{xi)ipi{x2)) = iZ{xi;x) X Gi{x,X2)- (65) 

This relation was obtained in Crocce & Scoccimarro 
(2006b) from a diagrammatic approach. Thus, we see that 
the three approaches (^5]) . ((55|) . and ([S^ are closely re- 
lated. In the integral method we simply absorb the matrix 
O into the response field x- Next, we can either take the 
limit rji — > —00 from the start, as for the action 5, or keep 
a finite rji in the computation, as for S. Then, we can take 
777 ^ — cxD in the final results for the nonlinear two-point 
correlation. Note, however, that for the approach of Crocce 
& Scoccimarro (2006a, b), which corresponds to the action 
S, it is not possible to take this limit in a practical manner, 
since one needs to keep track of the response TZ, which has 
no finite limit for r;/ — > —00. This leads to somewhat more 
complicated expressions than for the approaches based on 
the actions S and S of Eqs. (|i5|) . ([55|) . where the response 
functions R and TZ remain well-defined for r]j —00. Of 
course, the analysis described above also applies to the case 
of the gravitational dynamics. 

4. Large- expansions 

The path integrals ([45]) , ((54|) , and ((60)) can be computed by 
expanding over powers of the non-Gaussian part (i.e. over 
powers of Kg or Kg). This actually gives the usual pertur- 
bative expansion over powers of the linear power spectrum 
Pl (see also Valageas (2001, 2004) for the case of the Vlasov 
equation of motion). On the other hand, these path inte- 
grals may also be studied through large- expansions as 
in Valageas (2004). We focus below on the differential form 
(P5|) . but the formalism also applies to the integral forms 
([54)) and ()60p . Thus, one considers the generating functional 
Zn [7, h] defined by 

Zn[j, h]= J [dV^][dA] e^b-'>+"-^-^['^.^]l, (66) 

and one looks for an expansion over powers of 1/A'^, taking 
eventually N = 1 into the results. As discussed in Valageas 
(2004), the large- expansions may also be derived from a 
generalization of the dynamics to N fields This yields 
the same results once we deal with the long-wavelength 
divergences that constrain which subsets of diagrams need 
to be gathered. 

The interest of such large-iV expansions is to provide 
new systematic expansion schemes that may show improved 
convergence properties as compared with the standard per- 
turbation theory. Besides, it is clear from Eg. ([66)) that the 
symmetries of the system (e.g. invariance through trans- 
lations) are automatically conserved at any order. These 
methods have been applied to many fields of theoretical 
physics, such as quantum field theory (e.g. Zinn-Justin 



1989; Berges 2002), statistical physics (e.g. study of growing 
interfaces described by the Kardar-Parisi-Zhang equation, 
Doherty et al. 1994), and turbulence (e.g. Mou & Weichman 
1993). They are closely related at lowest order to the so- 
called "mode-coupling approximations" used for critical dy- 
namics, liquids, or glassy systems (Bouchaud et al. 1996), 
and to the "direct interaction approximation" used for tur- 
bulent flows (Kraichnan 1961). Therefore, it is natural to 
investigate their application to the cosmological gravita- 
tional dynamics described by Eqs.([T|)-([3]), which are similar 
to the Navier-Stokes equations. In some cases (e.g. Berges 
2002), it has been found that, whereas the simplest pertur- 
bative expansions give rise to secular terms (which grow as 
powers of time), the 2PI effective action method derived 
from such a large- iV method (discussed below in Sect. 14.2) ) 
could achieve a non-secular behavior and display relaxation 
processes. Of course, the actual behavior of such schemes 
depends on the specific problem. It has already been shown 
in Valageas (2007) that, for the case of the gravitational dy- 
namics in the expanding Universe, the large-iV expansions 
indeed show a qualitative improvement over standard per- 
turbation theory at one-loop order, as they display bounded 
oscillations (for the steepest-descent method) or decaying 
oscillations (for the 2PI effective action method) for the re- 
sponse functions instead of the secular terms encountered in 
the standard perturbative expansion (which gives increas- 
ingly large powers of time D{tY at higher orders). In this 
article, we investigate whether this good behavior extends 
to higher orders in the case of the Zeldovich dynamics. 

We discuss below both "linear schemes", such as 
the standard perturbation theory or the steepest-descent 
method of Sect. 14.11 which involve expansions over lin- 
ear two-point functions, and "nonlinear schemes" , such as 
the 2PI effective action method of Sect. 14.21 which involve 
expansions over nonlinear two-point functions themselves. 
By expanding different intermediate quantities or different 
equations (derived from the same equation of motion) , one 
obtains different methods that also correspond to different 
partial resummations. 

4.1. Steepest-descent method 

A first approach to handle the large- limit of Eq. ([S5)) 
is to use a steepest-descent method (also called a semi- 
classical or loopwise expansion in the case of usual quantum 
field theory with h = i/N). For auxiliary correlation and 
response functions Gq and Rq, this yields the equations 
(Valageas 2004) 

0{x,z).Go{z,y) = (67) 
0{x,z).Ro{z,y) = Soix-y) (68) 
Roix,z).0{z,y) = Soix-y), (69) 

whereas the actual correlation and response functions obey 

0{x,z).G{z,y) = i:{x,z).G{z,y) + n{x,z).R^{z,y) (70) 
0{x, z).R{z, y) = 6D{x^y) + S(a:, z).R{z, y) (71) 
R{x, z).0{z, y) = 6D{x-y)+ R{x, z).E(z, y) (72) 

where the self-energy terms S and 11 are given at one-loop 
order by 

T,{x, y)=4Ks{x; xi,X2)Ks{z; y, Z2)Ro{xi, z)Go{x2, Z2) (73) 
n(x, y)=2Ks{x; xi,X2)Ks{y; yi,y2)Ga{xi,yi)Go{x2,y2)-{7A) 
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Note that Eqs. (p7|) - ([7^ are exact and that the expansion 
over powers of 1/iV only enters the expression of the self- 
energy (|73|) - (l74| . Here we only kept the lowest-order terms 
(see Valageas 2004 for the next-order terms). We also took 
the limit rji — > —oo so that terms involving A/ vanish. The 
comparison of Eqs. ((67|) - ((68)) with Eas.([34|) shows that the 
auxiliary matrices Gq and Rq are actually equal to their 
linear counterparts: 



4.3. Role of self-energy terms 



Go — Gl, Rq — Rl- 



(75) 



Next, substituting Go and Rq into Eqs. ([751) - ([7i)l . we obtain 
the self-energies at one-loop order. First, we can note that 
S has the same form ([55)1 as the response R, whereas 11 
is symmetric and has the same form (I3ip as the two-point 
correlation G. Then, a simple calculation gives 

Eo(a;i,a:2) = -ujl0{rii - 1^2) So O^i - 'k2) 

1 j +^ loo 



where we define uji ~ uj{ki) as 



uj{k) = kay with 



in 



dwPUw). (77) 



Here is the variance of the one-dimensional displacement 
field Slo (or of the one-dimensional velocity dispersion up 
to a normalization factor). On the other hand, H is given 
at one-loop order by 

Ho(xi,X2) = 5i,(ki -f k2)e2"i+2''^Ho(fci), (78) 
with 



Ho(fc) = 2 J dkidk2(5i5(ki+k2-k)Pio(fcl)PLo(fe) 



7ri7r2 7r| 



and 

a(ki,k2) + a(k2,ki) 

TTl = , 712 = /3(ki,k2) 



(79) 



(80) 



Then, the response R and the correlation G can be obtained 
by integrating Eqs. ([7ni) - ([7T|) . 

4.2. The 2 PI effective action method 

As described in Valageas (2004), a second approach is to 
first introduce the double Legendre transform r[-0, G] of 
the functional W = \nZ (with respect to both the field i/j 
and its two-point correlation G) and next to apply the 1 /N 
expansion to F. This "2PI effective action" method yields 
the same equations ([70| -([72 |) . and the self-energy shows the 
same structure at one- loop order as ([75)) - ([7H) where Gq and 
Rq are replaced by G and R: 

E(a;, y)=4Ks{x; xi,X2)Ks{z; y, Z2)R{xi, z)G{x2, Z2) (81) 
H(a;, y) = 2Ks{x; xi,X2)Ks{y; yi,y2)G{xi,yi)G{x2,y2)- (82) 

Thus, the direct steepest-descent method yields a series of 
linear equations that can be solved directly, whereas the 2PI 
effective action method gives a system of nonlinear equa- 
tions (through the dependence on G and i? of E and H) that 
usually must be solved numerically by an iterative scheme. 
However, thanks to the Heaviside factors appearing in the 
response R and the self-energy E, these equations can be 
solved directly by integrating forward over the time 771. 



From Eq. fTT]) we can see that the self-energy S plays the 
role of a damping term. Indeed, Ea. ([7T|) has the form 
dR/dr]i = E.i? so that large negative values of S are as- 
sociated with a strong damping of the response function 
(Exact details are somewhat more intricate since Eq. ([7T|) is 
actually an integro-differential equation.) This agrees with 
Eg. ((76)) which shows that the one- loop self-energy Eq be- 
comes large and negative at high fc as Eq cx: — fc^. Thus, the 
self-energy E encodes the loss of memory associated with 
the nonlinear dynamics. 

On the other hand, we can see from Eq. lffD]) that the 
self-energy H is associated with the continuous production 
of power due to nonlinear mode couplings. Indeed, we can 
see from Eos. ([7n| - ([7^ that the correlation G can also be 
written in terms of the response R as 



(76) G{xi,x2) = Rx Gi X + R.U.R^, 



(83) 



and we let iji —00. The physical meaning of Eq. (|83|) is 
clear. The first term on the right hand side means that the 
fluctuations at the initial time 777 are merely transported 
forward in time through the response R. This is the only 
nonzero term in the linear regime (with R = Ri hence G — 
Gl)- The effect of the nonlinear dynamics is to modify the 
transport matrix R and to add a second term to the right 
hand side of Eq. ([55]) . The latter has the meaning of a source 
term that produces fluctuations with two-point correlation 
Il{ri[,r]2) the times ir]i,rf2) that are next transported 
forward to later times (771,772) by the matrices R{rji,ri[) 
and i?^ (77^, 7^2). 

5. Running with a high-Zc cutoff 

In a recent paper, Matarrese & Pietroni (2007a) introduce 
another approach to studying the gravitational dynamics 
within the hydrodynamical framework. It is also based on 
a path-integral formulation. Although they use the integral 
form of the equations of motion, as in Sect. 13.21 we briefly 
describe in this section how this method may be applied 
to the path integral (US)) . First, from Eq. ([^51) we define the 
generating functional Z[j, h] as 



Z[j,h] = I [d^][dA] 



J.4>+h.\-Sli>.\] 



(84) 



where we have introduced the external source h. This allows 
us to obtain the correlations of the response field A through 
derivatives with respect to h. Next, following Matarrese & 
Pietroni (2007a), we add a high-fc cutoff A to the linear 
power spectrum PLo{k) by changing the kernel A/, which 
appears in the action S of Eq. (|46|) as 



A, 



Aa = 6'(A- fci)A/(a;i,a;2). 



(85) 



Thus, the Heaviside factor 6{K — ki) removes the linear 
power at high wavenumbers fci > A, and we recover the 
full system in the limit A — s- cx). Then, the idea proposed 
in Matarrese & Pietroni (2007a) is to study the evolution 
of the system as a function of the cutoff A. Therefore, one 
first looks for equations that describe the dependence on 
A. Second, one derives some approximation for these equa- 
tions, for instance by a truncation of some expansion, and 
finally solves these approximate equations from A = up 
to A = CX3. First, the dependence on A may obviously be 
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described through the derivative of Z with respect to A, 
which reads 



OA 2 



dk(5c(A-fc) 



Shi{k,r]j)Shj{-k, rji) 



.(86) 



Next, introducing the generating functional W of the con- 
nected correlation functions, 



W ^\nZ, R{xi,X2) 



dR 



6j{xi)5h{x2) 



(87) 



=h=0 



we obtain from Eq. (|86p the evolution of the response R with 
the cutoff A as 

5'"^Plo(A) 



dA 



ixi,X2) 



E 



&w S]j{A — w) 



5hi{-w,rii)Shj{--w,r]i)Sj{xi)Sh{x2) ' 



Here we use the property ([17|) : (A) (A A) = 0. Next, to 
make some progress, one needs to obtain an expression for 
the fourth-derivative W^^'^ . Of course, in generic cases this 
quantity is not known exactly and one must introduce some 
approximations. The usual procedure is to write a diagram- 
matic expansion for W, using the path integral expression 
and to truncate at some finite order. For the cubic ac- 
tion pS)) . the lowest-order contribution is associated with 
the diagram of Fig. [l] which gives 



Fig. 1. The first diagram of the expansion of the fourth 
derivative ly*-^-* as in Eq. ([55)) . The big dots are the vertices 
Kg and the lines are the two-point functions R or G. 




= {lplX2h>^4)c 



Sji5h2Sh3Shi 

= 12RR{K,RKs)RR+ .. (89) 
Using expression ^U\i of the vertex Kg, this gives: 

^^^{k;f^,,7j2)=4e^'"PLo{A)Jd^v5D{A~w) di^jdrj' 
xRhi'^ [k; m,v)Ri'^i2 [k; 77', V2)Ri'^i3 {w; 77, Vi)Rt'^t4 ("w; v', Vi) 



xi?y (k - w; T], Tj')j^, i(w, k - w)7|.,, (k, -w 



(90) 



Following Matarrese & Pietroni (2007a,b), we note that at 
lowest-order we may replace the response functions on the 
right hand side in Eq. (|90p by the linear responses, which 
do not depend on A since they do not depend on the linear 
power spectrum, see Eq. ([5(Il) . Then, in the limit rji —00, 
only the linear growing modes of Ri'^^i^ and Ri'^i^ give a 
non-vanishing contribution: 



dR 



^=4Plo(A) J dwdoiA-w) j j dry' 6"+" 

X Rhiii'^ {m,v)RLz'^i2 iv', V2)Rlij iv,v') 

X7^.,(w,k-w)7;..(k,-w). (91) 



Using Eqs. (PT|l - (|^^ for the vertices 7^ we obtain: 



dR 4^,2^ r^^iDl-J 



-Rl- 



(92) 



Starting from the initial condition R{A = 0) = Rl, this 
yields at A = 00, 



R = Rl 



1 



(93) 



which agrees with the usual perturbative result at order Plo 
(see Eq. (|155l) below). The running with A of the two-point 
correlation G (whence of the nonlinear power spectrum) is 
then obtained by taking the derivative with respect to A of 
Eq. ([55)) and again using a loopwise expansion for the self- 
energy n. In practice, Matarrese & Pietroni (2007a, b) use 
some ansatz for G to derive a linear equation that can be 
solved up to A = cxD. Although they refer to this approach 
as a renormalization group method, we can note that it is 
somewhat different from the usual renormalization group 
techniques. Although considering the evolution with a cut- 
off A, one does not look for a fixed point of renormalization 
group equations that would govern the properties of the 
system in some large-scale limit. 

Matarrese & Pietroni (2007a, b) notice that, if we pro- 
mote the linear response Rl on the right hand side of 
Eq. to the nonlinear response R, we obtain the response 
(|133p with a Gaussian decay at high k as the solution of 
this linear equation. Note that all the previous steps apply 
identically to the case of the gravitational dynamics, where 
this procedure again leads to the response (|133p . In this 
case, this expression is no longer exact, but it does show the 
expected damping into the nonlinear regime. This remark 
suggests that this procedure may provide a very efficient ex- 
pansion scheme for the response function. However, this is 
somewhat artificial. Indeed, to derive Eq. ([5^ from Eq. ipO]) . 
one makes use of the properties of the linear response to 
simplify the right hand side, so that substituting back the 
nonlinear response R is somewhat ad-hoc (although it is 
correct at lowest-order it makes the procedure not system- 
atic). Moreover, it is clear that one can apply the same 
procedure to any other scheme that gives an equation of 
the form dR/da — F[Plo, R, G] where a can be any vari- 
able among {k, 771, 772, A, ..}. (For the large- iV expansions of 
Sect.m it would be 771.) Indeed, at lowest order one can al- 
ways simplify as a linear functional of Rl such that one 
obtains the exact response (|133p by substituting Rl ^ R, 
since the right hand side must be consistent at lowest-order 
with dR/da evaluated for the exact response (|133p . Thus, 
the latter satisfies the equation 



dR 
da 



dluRL 1 d 



dc 



2 da 



R, 



(94) 



which implies that at lowest order one can always write 
dR/da = F[Plo,R,G] as 



d_ _ dlnRL 
da da 



Rl 



.(95) 



where the dots stand for higher-order terms over {{Di — 
D2),oj}. For the specific case a = A, Eq. ipSj) leads back to 
Eq. iP^ . For the large- expansion schemes, Eq. ipS)) would 
be Eq. (|7ip with a 771 , the left hand side corresponding to 
O.R and the right hand side to S.i? at lowest order. Then, 
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substituting R to Rl on the right hand side of Eq. ([55)1 , one 
recovers Eq. and the nonhnear response l\l'S3\i with the 
Gaussian cutoff. This reasoning also apphes to the gravita- 
tional dynamics, to which one adds high-fc approximations 
so that the response (|133p (or a variant) still applies; see 
also Sect. |9] below. Therefore, recovering the response (|133p 
in this manner does not imply that the underlying expan- 
sion scheme is very efficient. However, Matarrese & Pietroni 
(2007b) argue that for the running with A it is possible to 
derive a stronger justification of this procedure which still 
applies at higher orders. Then, for cases where additional 
arguments can be obtained, such techniques based on the 
evolution of the system with respect to some parameter a 
may provide useful alternative expansion schemes. We dis- 
cuss this method, based on the dependence of the system 
on a high-fc cutoff A, in Sect. [T^ below within the frame- 
work of a simple systematic expansion, and we find that 
it actually gives similar results to the 2PI effective action 
method. 

6. Exact two-point functions 

For the Zeldovich dynamics, all quantities of interest can be 
computed exactly, since we know the solution l|10p of the 
equations of motion. This makes the Zeldovich dynamics 
an interesting test of approximation schemes, since we can 
compare their predictions with the exact results. As the 
equations of motion in the form (|lip - (jl2p are very similar 
to those associated with the exact gravitational dynamics, 
we can expect that the behavior of various approximation 
schemes will be similar for both dynamics. Therefore, we 
compute the exact two-point functions associated with the 
Zeldovich dynamics in this section. 



6.1. Two-point correlation 

As is well known, the two-point correlation G for the 
Zeldovich dynamics can be computed exactly from the so- 
lution pO| of the equations of motion (e.g. Schneider & 
Bartelmann 1995; Taylor & Hamilton 1996). Indeed, start- 
ing from the uniform density p at f — s- 0, the conservation 
of matter gives, before orbit-crossing. 



p(x)dx = pdq whence 1 -I- <5(x) 
This also reads from Eq. fTU)) as 



det 



9x 
dq 



= / dq(5n[x-q-s(q,77)] - 1, 



(96) 



(97) 



where s(q, 77) = £>+(r/)sio(q) is the displacement field. 
Note that this expression remains valid after shell cross- 
ing: all particles of Lagrangian coordinate q that happen 
to be at location x at the time of interest contribute to the 
right hand side. In Fourier space we obtain, for fc 7^ 0, 



dq 

(27r)3 



-ik.(q+s) 



(98) 



Therefore, the density-density two-point correlation reads 
as 



<5D(ki +k2)Gii(fci 
J (2^)6 



qi+k2.q2) 



i(ki.Si+k2.S2~ 



(99) 



Since the displacement field slo given by Eq. (l25|) is 
Gaussian, the average in Eq. (|M|) reads as 

^g-i(ki.si+k2.S2)^ _ ^-^{{kiiSii+k2iS2i)(kijSij+k2jS2j)) (100) 

where we sum over the 3D components i,j = 1, 2, 3. Let us 
define the displacement correlation ^^o: 

*o;»j(qi,q2) = (sL0ri(qi)sL0u(q2))- (loi) 

Thanks to statistical homogeneity, it obeys 

*o;»j(qi,q2) = *ory(qi - q2), (102) 

and from Ea. (|25p it is given by 



(103) 



vl^0;., (q)= J dke* '>^Pio(fc). 



Then, Eq.dM]) writes as 



Gil (fc; 771,772 



dq 



(2^)3 

^ ge''i + ''2fc,fe^.[*0..j(q)_cosh(,7i-»,2)*0;.j(0)]^ ^^Q4^ 

Using Eq. (|103p this can also be written as 



Gii(fc;77i,?72) 



dq 

(27r)3 



-zk.q 



X e'^''^^''^ ■f '^^ ''''t"' -PLo(»")[cos(w.q)-cosh(r)i--?72)] (105) 

The integration over angles in the exponent of Eq. (|105p can 
be performed analytically. Thus, let us define the quantity 
/(q;k) by 

I{q;k)^k,k,^o-.A<l)= J dwe'-^-^^^^^PLoiw). (106) 

Then, by expanding the exponential over spherical harmon- 
ics, we obtain 

I{q;k) = kH,{q) + e{l-3^l^)I2{q), (107) 



kq 



where we introduce 
47r 



In{q)^^l dw PLo{w)jn{wq), 

•J Jo 



(108) 



and jn is the spherical Bessel function of order n. Note 
that the variance ct^ of the one-dimensional displacement 
field, defined in Eq. ([77|) . also satisfies — /o(0). Then, 
Eq. (|105p reads (using the linear growth factor D ~ as 
the time-coordinate) as 



Guik;Di,D2 



X e 



dq 

(27r)3 

DiD2k^lIo + (l-3ti^)l2] 



cos(fc<7/i) 



(109) 



Following Schneider & Bartelmann (1995), we can perform 
the integration over angles by expanding the exponential 
and using the property 

r dfi cos(fcgM)(l - /i')" - nl (^A^ j„(fcg). (110) 



This gives 
Gii(fc;Z?i,Z?2)=e 



'^1 1^ pDiD2fe^(/o-2/2) 

27r2 



X > I UTU2—r— ) jn{kq). (Ill) 
n=0 ^ 
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6.2. Asymptotic behavior 



6.3. Response function 



From the explicit expression ()109|) . we can obtain the 
asymptotic behavior of the two-point correlation function 
in the highly nonlinear regime. Thus, we can formally write 
for a power-law linear power spectrum 



F[Alik;Di,D2)] (112) 



with 



F(:E)=y 0cos(qA*)e^/''""""'^"^""^"'+^'"'^'^^'^''^'"^l (113) 

where we make the change of variables q — )■ q/k,w kw, 
for 



Pwik) = 



whence 



1 f kV 



inkQ \ko 



Ai{k;Di,D2) = DiD2[^ 



n+3 



(114) 



(115) 



First, we note that infrared (IR) divergences appear in the 
one-dimensional velocity dispersion (defined in Eg. ([771) 1 
for n < —1 at low k, and in the integral over w in 
Ea. (|113p for n < —3. As is well-known, the IR divergence at 
n < — 1 should disappear for equal-time statistics (Vishniac 
1983; Jain & Bertschinger 1996) because of Galilean in- 
variance. This is explicitly checked in Eq. (|112p since for 
Di = D2 the prefactor of k^a^ vanishes so that the con- 
tribution associated with cry cancels out. Thus the equal- 
time nonlinear power spectrum is well-defined for n > — 3. 
Second, we can see that both and the integral over w 
in Eq. fTT^ diverge if n > -1 at high k. Thus, this UV 
divergence remains untamed in the full non-perturbative 
result (|112p . This is a qualitative difference with the true 
gravitational dynamics where such UV divergences are ex- 
pected to disappear in the exact nonlinear power spectrum 
for — 3 < rt < 1. However, this may require going beyond 
the single-stream approximation. Therefore, in the follow- 
ing we assume — 3 < n < — 1. After performing the integral 
over w and making a change of variable, we obtain 



F{x) = x"+^ — J dqq^ J d/icos ^a;"+ig^^ 



oc /•! 

dq q^ 
.^0 



xexp 



(-n-l)r[(4-n)/2]' 



(116) 



which shows that F{x) is well-defined for —3 < n < — 1 
and obeys the asymptotic behavior 



F{x) 



a;"+i for x^ \. 



(117) 



Thus, the equal-time power A^(fc;D) decreases in the 
highly nonlinear regime as 



A\k-D)^Al{k-D) 



for A? > 1. 



(118) 



Therefore, if Pwik) ^ fc" at high k the nonlinear power 
decreases as a power law P{k) ^ f^-3+3{n+3)/(n+i) ^j-^g 
highly nonlinear regime. 



Using the exact solution pOj) we can also compute the exact 
response function R. First, we note that, since the velocity 
field is decoupled from the density field, we have the simple 
exact result: 



i?21 = 0. 



(119) 



Of course, the linear solution (I30p is consistent with 
Ea. pi9p . Next, we can compute the response R12 of the 
density to a velocity perturbation as follows. At time ry^, 
before the velocity perturbation localized at time 772, the lo- 
cation and velocity of the particle of Lagrangian coordinate 
q are from Eq. pT ^ 



(120) 



X2 = q + L'2SLo(q), V2 =D2SLo(q), 



whereas at time 77^, after the velocity perturbation C2(x), 
we have: 



(121) 



= X2 , V+ = V2 - hT-(-2yJ-C2 



where we have used the definition of ■02 in Eq. p3l 
Therefore, the location of the particle at time ryi > 772 is 



xi = q + L>iSLo(q) 



D2 



(122) 



The density contrast is again given by expressions of the 
form (f96 |) - ([98|) so that ■!/'i(ki, 771) = (5(ki,77i) reads for fci 7^ 
as 



■0i(ki,77i) 
Definition 
-Ri2(ki,77i;k2,772 



(27r)3 

of the response function reads here as 
<5V'i(ki) 



^ <5C2(k2) 

Then, using the expression 
-V-i.C2 = j dke'''-^7^,C2(k) 
we obtain from Eq. p23p 



(124) 



(125) 



R12 



Di-D2 ki.ka 



dq 

D2 kl J (27r)3 

X ^g-«(-Diki--D2k2).SLo(q)^ 



„i(k2-ki).q 



(126) 



Because of homogeneity, the average in Eq. (|126p does not 
depend on q so that the integral over q yields a Dirac factor 
(5D(ki — k2). On the other hand, since s^o is Gaussian the 
average can be easily performed as in Sect. 16. 11 which gives 
for Di> D2: 

R,2ik; D„D2) = e-'^iD.^^^fk^.l (127) 



i?Li2e 2 



(128) 



where Rli2 is the linear response from Eq. (j30p and we have 
introduced w(fc) = ka^ as in Eq. fTT]) . The computation of 
Rii proceeds along the same lines. A perturbation Ci(x2) 
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of the density field at time 772 does not modify the velocity 
field, and we obtain, for ki ^ 0, 



V'llki,??!) 



dX2 



dx2 



[1 + Si^2,vt)] 



-iki .xi 



(27r)3 



-zki .xi 



dq^Z3[x2-q--D2SLo(q)] +6(^2) 



(129) 



where we have used Eq.([97|) and xi,X2 are the locations 
at times rji , 772 of the particle of Lagrangian coordinate q. 
This gives 



Rii = (det 



/9X2 

\dq 



-j(Di-_D2)k.Sio(q) 



Expanding the determinant, we find that most terms cancel 
out, and we obtain the simple result: 



In a similar fashion, for i?22 we can write: 



(131) 



M^i,m)-I ^det 



-iki .xi 



dk 



ki.k 



(132) 



The computation is slightly more intricate than for Rn, 
since Xi also depends on C2; however, most terms cancel out 
and we recover again the same form as in Eas. (|128p . (|13ip . 
Thus, the exact nonlinear response function is merely given 

by 



i?(fc;i^i,i52)-i?Le-5(^i-^^)^ 



(133) 



that is, all linear components are multiplied by the same 
damping factor. We can see from Eq. p33p that the re- 
sponse function only depends on the linear power spec- 
trum through the linear velocity dispersion a^, and on scale 
through = k'^<Jy- This property extends to the self- 
energy E which is related to R through Ea. ([7T|) . This is a 
big simplification with respect to the gravitational dynam- 
ics, where the response R and the self-energy E depend on 
the detailed shape of PLo{k), see Valageas (2007). However, 
even in that case it appears that the behavior of the re- 
sponse function is mostly governed by = fc^cr^; see for 
instance the analysis in Sect. 5.2 of Valageas (2007). This 
also shows that both dynamics share important features. 

6.4. Damping self-energy S 

The self-energy E introduced in Sect. |4] is usually obtained 
as a series of diagrams from the path integral (j66p . However, 
since we know the exact response function R, we can di- 
rectly compute E from Eq. (|7ip . From the simple result 
(|133p . we can see that the matrix structure of R, hence of 
E, is not changed by the nonlinear corrections. Therefore, 
from Eq. ((76)l and Ea. (|133[) . we write the exact self-energy 
E as 



^{k;Di,D2)^^o<jHDi-D2)l 



(134) 



where the matrix Eq was obtained in Eq. (|76p . To generalize 
the calculation for future use, we consider a response of the 
form 

R{k;Di,D2)^ RLr{t), t = uj{Di-D2), r(0) = 1, (135) 



where the constraint r(0) = 1 comes from Eq. 
Substituting Eqs. (fTM|) - P^ into Eq.dTTl) yields 

r'{t) = - [ dt' a{t-t')r{t'). 



(136) 



The fact that the system ((7T|) can be reduced to Eq. (|136p 
shows that the scalings (I134p - (|135p are indeed self- 
consistent. Note that the functions r{t) and a{t) are defined 
for i > 0, because of the Heaviside factors 9{rji ~ 772) in R 
and E. Then, introducing the Laplace transform 



(130) ^(s) 



At e-''*r(t). 



we obtain from Eq. (|136p 
sf(s) — 1 = —a[s)r(s). 

For the exact nonlinear response (|133p . we have: 



T(t) 



-t'/2 



, f(s) 



which gives 
a{s) 



7rerfc(s/\/2) 

where erfc(x) is the complementary error function: 



erfc(a;) — —= 



dt e-* 



(137) 



(138) 



(139) 



(140) 



(141) 





Fig. 2. Left panel: the functions r{t) and a{t). Right panel: 
the function a{t) multiplied by a factor e* 1"^ . 



On the other hand, if we write the expansion of a{t") 
around i = as 



t > : G(t) = ^(-1)V, 

p=0 



{2p)V 

we obtain by substituting into Eq. (|136p 
p 

dp = (2p + 1)!! - ^P-rn{2m ~ 1)!!, 

m— 1 

and the first few coefficients are 

do = 1, (Ji = 2, (T2 = 10, (73 = 74, (T4 = 706, ... 



(142) 



(143) 



(144) 
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We note that this series also appears in other problems of 
field theory as the number of Feynman diagrams associ- 
ated for instance to a cubic complex action with two fields 
(Cvitanovic et al. 1978). This is not surprising since in our 
case we also have a cubic action P5|) with two fields tp,X. 
We show in Fig. [2] the behavior of the self-energy function 
a{t) computed from the expansion (I142|) . We can see that 
it shows a fast decay, together with oscillations (but the 
numerical range is too small to check whether the asymp- 
totic is of the form e~* cos(t)). Thus, the self-energy ex- 
hibits a more intricate behavior than the response r{t). 
This may explain why path-integral methods based on the 
Schwinger-Dyson equation (|7ip have difficulty reproducing 
the response R from simple approximations to E, as we 
shall see below. 

6.5. Self-energy II 

We can also compute the self-energy IT from Eq. lfSS]) . which 
gives 



n(a;i,a;2) = (O - I]).G.(0 - E)^ 



(145) 



using Eas. ([7T|) - ([7^ . Then, using the exact expressions of 
the two-point correlation G and of the self-energy S, we 
can obtain 11. However, using the exact two-point correla- 
tion Gij would give intricate expressions, as the velocity 
correlations involve a few prefactors up to order P|q in 
front of the exponential of Eq. (|105p . Here we are mostly 
interested in the qualitative behavior of the self-energy H, 
therefore, we use the simple approximation 



G~G with G = Gil 



1 1 
1 1 



(146) 



where Gn is the exact density-density correlation derived 
in Sect. 16.11 The simple form (|146p is also consistent with 
the linear regime limit (|27p . Then, expanding the exponen- 
tial of Eq. (|109p . we can write 



G(fc;Di,D2) 



dq 

(27r)3 



n=l 



I{q-kY 



xG„(fc,Z?i)G„(fc,Z?2)' 



(147) 



where /(q;k) is defined in Eq. (|107p . and we introduce the 
vectors G„ defined by 



Gnik,D)=e- 



(148) 



In Ea. (|147p we use the fact that the term n = does not 
contribute to G{k). Then, the self-energy H associated with 
G through Ea. (|145p reads as 



n{k;Di,D2) 



with 



dq 

(27r)3 



1 



n=l 



/(q;k)^ 



xn„(fc,z?i)n„(fc,i?2) 



H„(fc, D)^(P- E).G„ = ^„(fc, D) \ 



(149) 



(150) 



Using Eq. (|134p we obtain 



'■D^ f dD'a[uj{D - D')]D''^-^e 
Jo 



(151) 



Then, using Eq. (|136p we can check that tti = 0. Moreover, 
if we apply the response i? to H as in Eg. ([55)1 . we obtain 
obviously from Eq. fTSOll and Eq. ([72|) 



xG„{k,Di)Gnik,D2f, 



(152) 



and we recover the full two-point correlation G by noticing 
that the term n = 1 in Eq. (|147p is equal to i? x G/ x 
(with rji —f — oo) in agreement with Eq. 

7. Standard perturbative expansions 




0.01 



0.1 1 

k [h Mpc-i] 



Fig. 3. The standard perturbative expansion of the re- 
sponse function over powers of Plq as in Ea. (|155p . We only 
show the density-density component i?ii for clarity. The 
solid curve i?NL is the exact response (|133p , whereas curves 
labeled R'^p^ are the expansion of the response function up 
to order p over Plq. Here we consider a ACDM Universe 
with "n = —2" normalized as in Ea. (|153p . but the results 
are identical for any CDM cosmology up to a rescaling of 
k. 

In this section, we describe the usual perturbative ex- 
pansion over powers of the linear power spectrum Plq ap- 
plied to the Zeldovich dynamics. In this article we are con- 
sidering a ACDM universe with f^m = 0.3, 17a = 0.7, ilb — 
0.046,0-8 = 0.9, a reduced Hubble constant h = 0.7, and 
we use the linear power spectrum given by the CAMB 
Boltzmann code (Lewis et al. 2000). This gives at z = 
for a smooth linear power spectrum taken from Eisenstein 
& Hu (1998): 

^Lo(^o) = l, t^(A:o) = l-3, for /sq = 0.21/iMpc"\ (153) 
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However, until Sect. 1 131 where we focus on weakly nonlinear 
scales and baryonic acoustic oscillations, we use a power- 
law linear power spectrum (|114|) with n = —2 (except in 
Sect. [TU| . normalized as in Eq. (|153p : 



Aio(fc) = k/ko, 



(154) 



R^P\k;D,,D2)^RLj2 



i-iy 



ml 



,(155) 



0.1 



- 1 1 1 1 1 1 1 1 II 
Zj = Zp = /' 




AZ(i)/' 




/ </\ ' \ 1' 
//.-■' ^ ( 




ji 

1 I I '! 


; A2<2> \ 


0.1 


1 



k [h Mpc-i] 



trative purposes we consider the power-law linear power 
spectrum (jl54p . In this case we can perform the integrals 
in Eq. (|116p and for the equal-time nonlinear power we ob- 
tain: 



where fco is given in Eq. (|153p . This is not important for the 
response function R, which only depends on Lo{k) — ka^, 
whatever the linear power spectrum, but this will allow us 
to simplify the analysis for the nonlinear two-point cor- 
relation G. First, we consider the response function R. 
Expanding Eq. (|133p over powers of Plo is equivalent to 
expanding over powers of w^, since oc Plq from Eq. ([77)) . 
Therefore, the response function R^p^ expanded up to order 
Plo is: 



A2 = F(A|) with F{x) = xF{xtt/8) 
and 

3y(l + yTTp')Arctan 

Hy) 



(156) 



This expansion converges absolutely at all times and on 
all scales, but the convergence is not uniform. Thus, the 
convergence rate is very slow for uj{Di — D2) 1 and for 
any finite order i?^^') grows without bound at large times or 
wavenumbers instead of decreasing. From Eq. (I155|) we can 
see that in order to obtain a reliable prediction at a given 
scale we need to go at least up to order p ^ u}'^{Di — Dq)^ ■ 
We display the first few terms in Fig. [31 which clearly shows 
that increasing the order p improves the agreement with 
the exact result on weakly nonlinear scales but worsens the 
prediction in the highly nonlinear regime. 



Fig. 4. The standard perturbative expansion of the two- 
point correlation over powers of Plo from Eq. (|158p . We 
only show the density-density equal-time logarithmic power 
A2(fc;_D) at redshift z = 0, for the case n = —2. The solid 
curve A^L is the exact power from Eqs. (|156p - (|157p . whereas 
curves labeled A^^^'^ are the expansion up to order p-\- 1 over 
Plo- Higher-order terms grow as fc^"*"^ at high wavenumbers. 
The perturbative expansion diverges beyond the vertical 
dotted line, for fc > 0.53/iMpc"^ 



Next, we turn to the standard perturbative expansion 
of the two-point density-density correlation Gn. For illus- 



^"^^y^y 4(1 + ^2)5/2^^2 + + 2VTTy2 

3y(-l + + y2)Arctan 



+ - 



4(1 -I- y2)5/2y^2 + y2 _ 2VTTy^ 



(157) 



Note that the last two terms have not been correctly written 
in Taylor & Hamilton (1996). The expansion over powers 
of Plo corresponds to the expansion of F{x) over powers of 
x, and we obtain 



A2.Ai + ^Ai-^Ai-i^. 
^ 64 ^ 32 ^ 8192 



(158) 



Contrary to the response function, we can see from Eq. (jl57p 
that this expansion diverges for A|^ > S/tt. Therefore, 
one cannot describe nonlinear scales from this perturbative 
expansion, and going to higher orders only improves the 
predictions for weakly nonlinear scales where A| < 8/7r. 
We can see from Eqs. (jll3p - (jll6l) that the radius of con- 
vergence of the perturbative series is zero for n < —2. 
(Since at large q we encounter an integrand of the form 
/ dge~^' " , which gives rise to a singularity at a; < 0.) 
For a ACDM cosmology, the slope of the linear power spec- 
trum goes to 71 = 1 on large scales; therefore the pertur- 
bative series should always converge. However, on small 
scales where n < — 2 the perturbative expansion is likely 
to be useless (except for the quasi-linear regime) since the 
series only converges because of the behavior of the lin- 
ear power spectrum on unrelated large scales. We compare 
the first few terms A2(p) of this perturbative expansion 
with the exact nonlinear power in Fig. 01 In agreement 
with the analysis above, we can check that the perturba- 
tive predictions provide a good match on quasi-linear scales 
A| <C 8/7r, fc ^ 0.53/iMpc^^ and becomes useless deeper 
into the nonlinear regime. 

Crocce & Scoccimarro (2006a) notice that instead of the 
standard perturbative expansion (|158p . one could keep the 
exponential factor e~^^^~^^^^ " ^2 jn expression (|109p and 
only expand the last exponential e^^^^^ . In this fashion, 
all terms of the new series are positive and damped by the 
exponential factor at small scales, which gives a seemingly 
better behaved expansion. We display the results we obtain 
for the case n = —2 in Fig. [51 when we consider the series 



A2 = 



37r2 



64 



(159) 



where we use the normalization of Eq. p53p for uj{k). Of 
course, for a power-law linear power spectrum, we actually 
have u — 00, but Eq. (|159p describes a CDM-like power 
spectrum, such that n = —2 on the weakly nonlinear scales 
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Fig. 5. The perturbative expansion of the two-point cor- 
relation over powers of Plo from Eq. (|159p where we keep 



apart. 



The sohd curve A^^ 



the exponential factor e 
is the exact power from Eqs. (ll56p - (ll57p . whereas cm^ves la- 
beled A^'^P) are the expansion up to order p +\ over Plq- 
All terms decay as "-"^ at high k. The perturbative 
expansion again diverges beyond the vertical dotted line, 
for k > 0.53/iMpc-^ 



of interest and w is made finite and normalized to Eq. pSSp 
through IR and UV cutoffs. Fig. [5] shows that indeed the 
terms at each order over Plo are positive, and the series 
looks better behaved. However, it is clear that the series still 
diverges beyond k ~ 0.53/iMpc~-^ as for (|158p . Moreover, 
we note that increasing w(fco) (by moving the IR and UV 
cutoffs) would move the Gaussian cutoff towards smaller k 
and would worsen the agreement with the exact nonlinear 
power. In fact, for general power-spectra where the linear 
velocity dispersion is not necessarily governed by the 
weakly nonlinear scales of interest, the rewriting associated 
with Eq. (jl59p does not always improve the agreement with 
the exact power. Power-law linear power spectra are a clear 
example of such cases, since w = c» so that the expansion 
(|159|) is not well-defined (unless one absorbs IR divergences 
into a renormalized crJ3) , even though the equal-time power 
remains finite, and the usual perturbative expansion (|158p 
is well-defined. Note that the expansion schemes based on 
Ea. l|83p . through path-integral or diagrammatic methods, 
do not exactly correspond to the expansion (|159p . Although 
the first term Rx Gj x may exhibit the factor , 
this is no longer true for the second term (because of the 
integrations over time in R.Il.R^), as shown by the discus- 
sion in Sect. 18.31 We note, however, that since the response 
function R also involves the quantity lo we can expect such 



^ The procedure (|159|) could be modified to apply to cases such 
as n = — 2 without any IR cutoff, as the exact expression (|112p 
shows that the IR divergences are fully absorbed into a; — ka^. 
Then, by splitting all IR-divergent integrals into a cr„ part and 
a finite part and using a given value for (t„, one obtains finite 
expansions. The value of is irrelevant for exact equal-time 
statistics (but would remain in approximate quantities obtained 
as in Ea. (|159p by truncation at some order). 



schemes to fail in cases where the velocity dispersion cr^ is 
governed by scales that are very far from those of interest. 



8. Direct steepest-descent method 

8.1. Response R and damping self-energy E 




0.1 1 

k [h Mpc-i] 



Fig. 6. The expansion of the response function defined by 
the direct steepest-descent method, from Ea. (|I62p . We only 
show the density-density component Rn for clarity. The 
solid curve i?NL is the exact response (|133p . whereas curves 
labeled R'-p^ are the expansion of the response function up 
to order p. 



We now investigate the properties of the steepest- 
descent method described in Sect. 14. ll We first consider the 
damping self-energy E and the response R. The self-energy 
E can be obtained at one- loop order from Eg. ([75)1 . which 
gives Eq. (|76p . and at higher orders by including higher- 
order diagrams. However, for the Zeldovich dynamics, it 
can be directly obtained from the exact expression derived 
in Sect. 16.41 From Ea. (|134p it is clear that series (|142p cor- 
responds to the perturbative expansion of the self-energy E 
over powers of Plq (through powers of to^). The one-loop 
result ((76|) simply corresponds to the first term cr^^^(i) — ctq 
(whereas the linear regime corresponds to a — 0) because 
of the prefactor in Eq defined in Eq. (|75|) . In this fashion 
we obtain the self-energy E up to order p as 



E(rf = Eoa(rf = Eo g(-l)'"a„.M^i-:^. (1 



m— 



60) 



This in turn determines the response function through the 
steepest-descent method described in Sect. 14. 1] as 



R^PHk;Di,D2)^RLr<^''^HDi~D2)], 



(161) 



where AP^t) is obtained from a^P\t) through Eq. ()136p . 
From Eq. (fT5E)) this gives 



o2p-l 



+ Er=o(-i)'"^™*''^"'""^ 



(162) 
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For the first few terms this gives 



1 2 

j.{o) ^ 2, r^^) = cost, r(^) = -cosht + - cos \/2t. 

o o 



(163) 



Indeed, the perturbative expansion of the self-energy E 
gives a(t) as a series over powers of t, whence a{s) as a 
series over powers of 1/s. This yields a rational function 
for r^P\s) and a sum of exponentials for r'-P\t). At or- 
der p — I, the arguments of the exponentials are imagi- 
nary (they are given by the roots of the denominator of 
f^P\s)), but real parts appear at order p = 2, 3 (and pre- 
sumably at higher orders). Since the denominator is even, 
see Eq. (|162p . the roots appear by pairs ztsi so that both de- 
caying and growing exponentials appear. Therefore, the di- 
rect steepest-descent method cannot reproduce the decay of 
the response function at large t (i.e. in the highly nonlinear 
regime). Of course, at a given order p, the response A^^t) 
agrees with the expansion of the exact response (jl55p up to 
order t'^P. We display the first few terms i?^^-' in Fig. [HI The 
behavior is actually rather close to the standard perturba- 
tive expansion shown in Fig. [31 as the higher orders slowly 
improve the agreements over weakly nonlinear scales but 
explode increasingly fast into the highly nonlinear regime 
(but their amplitude grows even faster as exponentials in- 
stead of power laws). 



8.2. Pads approximants 
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Fig. 7. The expansion of the response function defined by 
the Fade approximants from Eqs. (|17ip . At each order, one 
obtains a sum of cosines, with an offset for even p. 



The remarks above suggest that we may improve the 
expansion (|162p by looking for a Fade approximant to the 
rational function f^^^ (s), which has the same expansion over 
1/s up to 1/s^P"'"^. In fact, since we know the exact response 
function, we can directly obtain the series of Fade approx- 
imants from Eq. (jl39p . First, from the expansion of e~* 



at i = 0, we obtain formally the expansion of f(s) over 1/s 
as 



p=0 



(164) 



Note that this only provides an asymptotic series for f(s) in 
the limit s — )■ oo. Here it is convenient to make the change 
of variable y = 2/s^ and to define f{y) by 



1 



f(s) = -r{y = 2/s ), r{y) 



OO 

E 

p=0 



{-Ifrpy". 



From Eq. (|164p we obtain for the coefficients rp 



(2p-l)!! rb+1/2] 



2P 



dt 



TTt 



tP. 



(165) 



(166) 



This shows that the expansion (|165p is a Stieltjes series 
since the coefficients fp are the moments of a real positive 
function defined over t >0 (here of the function e^YV^rt), 
see Bender & Orszag (1978). Since Carleman's condition 

is fulfilled, J2^p^^^^ — function r{y) is uniquely 

determined by its asymptotic expansion (|165p . which can 
be resummed as 



riy) 



dt 



l + yt ^f¥t 



(167) 



Then, all coefficients of the continued-fraction representa- 
tion of f{y) are nonnegative, and both Fade sequences 
and Ppj^x converge monotonically to f(y) as 



Pl{y) < pKv) <Pi{y)<-< r{v) 
r{y) < ■■ < Piiy) < Pliv) < Poiy) 

and 

f{y) = hm P;+i(y) - lim PP{y). 



p — >-oo 



p — !-00 



(168) 
(169) 

(170) 



By contrast, the usual perturbative expansion (|155p . which 
is associated with expansion (jl64p . amounts to approxi- 
mate r(y) by a polynomial, that is, by the sequence Pq (y), 
whereas the steepest-descent method (|162p amounts to ap- 
proximate l/r(y) by a polynomial, that is, f{y) by the se- 
quence Pp{y)- As we have seen above, these two sequences 

do not converge very well since they give a response r^^^ (t) 
that grows without bound for i — > oo. On the other hand, 
the sequence associated with (|170p gives 



f^°\y) = Po° = l, fWiy) = P? = 



1 



l + y/2' 



1 



l + 3y/2' 



(171) 



The first two terms give the same results r*^°^ and A^'' as 
Eq. (|163p but the next two terms give: 



,(2) 



it) 



r-(3)(t) = 



2 

3 + 
V6- 



2^6 

^/6- 



cos v 3 — V6t 



+ - 



2VG 



Vet, 



(172) 



(173) 
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which do not grow exponentially at large t any more. 
Because both the sequences and Pp+i and f(y) are 
Stieltjes functions, they are analytic in the cut plane 
I arg(y)| < tt, and all poles of the Pade approximants PP 
and -Pp+i lie on the negative real axis (Bender & Orszag 
1978). Therefore, Ea. (|165p shows that all poles of the as- 
sociated rational function f^P\s) lie on the imaginary axis 
and appear by pairs ±isj together with a pole at s = 
for even p. Then, the response factor A^^t) is a sum of 
cosines cos(sji) plus a constant for even p, in agreement 
with Eqs. (|172p - (|173p . This is a clear improvement over the 
standard expansion (|155p and the direct steepest-descent 
expansion p62p - (|163p . However, even this expansion can- 
not recover the Gaussian decay e~* which is replaced 
by fast oscillations. Nevertheless, this gives rise to an effec- 
tive damping (for odd p) once the response function is inte- 
grated with some weight function. We compare the first few 
terms to the exact nonlinear response in Fig. [71 Of course, 
we recover the same agreement at low k as with the other 
expansions displayed in Figs. [3] and [SI but as explained 
above, the response remains bounded at high k with fast 
oscillations. 



8.3. Correlation G and self-energy 11 

We now consider the predictions of the direct steepest- 
descent method for the two-point correlation G and self- 
energy n. As in Sect. [7l we consider the case of a power- 
law linear power spectrum n = —2, which simplifies the 
calculations, since integrals over wavenumber k have al- 
ready been performed following Eqs. (|156p - (|157p . Moreover, 
as in Eq. (|159p and Fig. [51 we consider a finite linear ve- 
locity dispersion ijj{k) = ka^ normalized in Ea. p53p . As in 
Sect. 16.51 in order to simplify the computations, we approx- 
imate the matrix form of the correlation Gij by Eq. (|146p . 
Thus, we keep the exact nonlinear, density-density correla- 
tion Gil, and make approximation Gtj = Gn for all 

Using Eq. (|112p and expanding the exponential e^^^'^^ and 
F{x), we write Gn as 



F, 



47rfc3 



p\m) 



where Fp is the p-th derivative F{x) dX x ~Q: 



p=0 



p\ 



(174) 



(175) 



The first few coefficients Fp are given in Eq. p58p . Note, 
however, that for the case n = —2, the Taylor series (|175p 
diverges for a; > S/tt as seen in Sect. [7l Then, we can write 
the matrix two-point correlation G as 



G 



F, 



^^A%u;^"'Gp+,n{Di)Gp+.n{D2f, (176) 



where the vectors G„ have been defined in Eq. (|148p . Thus, 
as compared with Eq. (|147p . the use of a power-law linear 
power spectrum has simply replaced the integral over q by 
a discrete sum. Then Eqs. (|149p - (|152l) still apply, once we 
replace the integral over q by this discrete sum. We now 
need to compute the functions 7r„(Z?) defined in Eq. (|15ip . 



Using the expansion p42p . we can perform the integrations 
over D' and obtain 



n-1 



i-iy 



m=0 

m!2'"+i " 



m + 1 



(n + 2m) 



m\ \ 2 

{n + 2i-iy. 



(177) 



Of course, from Eq. (|143p we can check that tti = 0, in 
agreement with the result already obtained in Sect. 16.51 
Substituting Eq. l|177p into the analog of Eq. (|149p . which 
reads here (as Eq. ()176p for G) as 

^=J^J2 ^Ai>2'"fl,+,„(i?i)IT,+„,(i52)^, (178) 



we obtain the self-energy 11 up to order P£q as 
u(p\k;D,,D2)^—^Y.^n{k;D,,D2)(l }), (179) 



71=1 



where n„ cx P^q^ and the first few terms are 



Hi = Ai ( ^Al+Luiuj2] 



U2 = -Al{Ljf+Ujl)(^Al+UJiUJ2] 



3_7r2 

64 y 

3_7r2 
64 
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where we have introduced (we recall that DiZ32A|q) 



(180) 



+AI{-^AI + ^AIloiu;2 + 2loIluI] (181) 



LUi — DiLu{k) = Dikffy, LU2 = D2i^{k). 



(182) 



As for the standard perturbative expansion ()158p applied 
to the correlation function G, the series expansion (|179p 
for the self-energy 11 gives higher-order terms that grow in- 
creasingly fast into the highly nonlinear regime. Moreover, 
we can expect that it diverges at equal times for A\ > S/tt 
as for G. 

Then, the direct-steepest descent method prediction at 
order p is obtained by applying Eq. ([551) using the response 
R^P^ and the self-energy 11^^) at that order. Note that this 
expansion does not have the form of a standard perturba- 
tive expansion such as Eqs. (|158p or (|159p . since all terms in 
the series get modified as we go to higher orders, because 
the response i?*-^-* has a different functional form, see (|163p . 
At hnear order p = we have 11 = and A^'"^ = A|. At 
first order p = 1, we obtain from Eq. (|163p 



A'^'\k;D,,D2 
'37r2 



A2 



64 



A| + CJitJ2 



^ cos(wi) cos(a;2) 

sin(cL;i) sin(cL;2) 



^1^2 



(183) 



Of course, if we expand Eq. (ll83p over powers of w, we re- 
cover the usual perturbative result, and at equal time we re- 
cover Eq. (|158p up to order A\. Moreover, we already know 
that at equal times Di = D2, all terms uj must cancel out 
up to order P£q ^ because the exact power (|156p does not 
depend on ui. This is a general result that does not de- 
pend on the shape of the linear power spectrum and can 
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k [h Mpc-i] 

Fig. 8. The expansion of the two-point correlation de- 
fined by the direct steepest-descent method, from Eq. (|179p . 
Higher-order terms display an exponential growth at high 
k for p> 2. 




0.1 1 



k [h Mpc-i] 

Fig. 9. The expansion of the two-point correlation defined 
by the direct steepest-descent method, from Eq. p79p . but 
using the Pade approximants (|172p - (|173p for the response 
function. Higher-order terms display a power-law growth at 
high k. 



be seen from Ea. (|105p . This is related to the cancellation of 
IR divergences recalled in Sect . 16 . 2l associated with Galilean 
invariance. In a similar fashion, at order p — 2, we obtain 
from Eq. p63)) 

+Al (-^Ai + ^^Aic.ic.,-H2c.M) rr'rf^^'' (184) 




k [h Mpc-i] 

Fig. 10. The expansion of the two-point correlation defined 
by the direct steepest-descent method, from Eq. (|179p . but 
using the exact response function p33p . Higher-order terms 
display a power-law growth at high k. 

where we have defined > 1) 

r[P^ =Ap\cj,), r'f''^ ^ f Att'~^ M[uj,{l-t)]. (185) 

We display our results for the first few terms in Fig. [51 
We can see that the exponential terms associated with the 
response function (see Eq. ()163|) ) make the higher-order ap- 
proximations grow exponentially at high k (for p>2). This 
very strong growth makes the series of little practical value 
for p > 2, where the convergence is not better than for 
the standard perturbative expansion displayed in Fig. 2) 
We note that the equal-time power can become negative at 
high fc, whereas both H and G should be positive matrices 
from Eg. ([55)1 . which implies that the components Gu and 
Tin are positive at equal times. This failure comes from the 
truncation of the self-energy H in Eq. (|179p . which breaks 
the positivity of H. Indeed, from Eq. (|18ip and the scalings 
A| c>c k and a; cx fc, we see that II2 ~ — 37r^/32A|^ijj^ at 
high k. 

We show in Fig. [5] the results obtained from the expan- 
sion ()179p for the self-energy H but using the Pade approx- 
imants ()172p - (|173p for the response function. The higher- 
order terms no longer grow as exponentials at high k but as 
power laws, since the response functions only contain con- 
stants and cosines instead of exponentials. However, this is 
not sufficient for significantly improving the convergence of 
the series to the exact nonlinear power A^. 

Finally, we show in Fig.[Tn]the results obtained from the 
expansion (|179p for the self-energy H when we use the exact 
response function p33p . Note that in this case the lowest- 
order approximation A^^''^ is equal to the linear power mul- 
tiplied by a Gaussian damping factor. It is actually equal to 
the first term of Eq. (|159p . Higher-order terms do not exhibit 
such a Gaussian damping because of the time-integrals in- 
volved in the last term of Eq. ([551) ■ Indeed, expansion (|179p 
yields power laws over Di,D2, which gives enough power 
generated at all times to build large density fluctuations 
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into the nonlinear regime. We can see from Fig. [TO] that 
using the exact response function is not enough to signif- 
icantly improve the convergence of the series obtained for 
the matter power spectrum. Therefore, it appears that to 
improve the results one should use other expansion schemes 
for the self-energy 11: improving the response function alone 
does not help much. 

8.4. Using a Gaussian decay for the self-energy U 





Zj = Z3=0 
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Fig. 11. The expansion of the two-point correlation from 
Ea. (|186p . using the exact response function (|133p and the 
expansion (|186p with Gaussian factors for the self-energy 
n. All terms display a Gaussian decay at high k. 

The previous discussion and the results obtained from 
expansion (|159p . shown in Fig. O suggest that it may be 
useful to factor out a Gaussian term of the form e^^ " 
from the self-energy 11. Therefore, we now replace expan- 
sion (fT79l) by 



n(fc;Z?i,i?2) = 



1 1 
1 1 



(186) 



where n„ is again of order P^q^ and uJi — Diuo as in 
Eq. (|182p . This expansion can be obtained as in Sect. 18.31 
or it can be derived from the expansion p79p by mul- 
tiplying by a factor e*^"i"'""2)/2 and again expanding over 
powers of Plq ■ We display in Fig. [11] the results obtained 
from Eq. (|186p using the exact response function (|133p that 
also shows a Gaussian decay at high k. Of course, because 
of the Gaussian damping factor introduced in Eq. (|186p all 

,2 2 

terms decay as e" at high k and the expansion looks 
better behaved than the one displayed in Fig. [TO] Besides, 
the approximation obtained at a given order seems to pro- 
vide good accuracy over a slightly wider range than in both 
Figs.[TOland[5l Thus, decomposing the two-point correlation 
in terms of response function and self-energy as in Eq. (|83p 
and using a Gaussian decay ansatz at high k appears to 
be a good scheme. However, this is somewhat artificial (see 
Sect.[l]). 



0.1 




k [h Mpc 



Fig. 12. The expansion of the two-point correlation from 
Eq. (|187p . using the exact response function (|133p and the 
expansion (|187p with Gaussian factors for the self-energy 

n. This Gaussian factor only damps the contributions at 
different times, and high-order terms display a power-law 
growth at high k. 

As explained in Sect. [9] below, the Gaussian damping 
factors ^ '^"/^ merely correspond to the linear displace- 
ment field that moves the location of large-scale structures 
between different times. However, this process does not af- 
fect the matter clustering, and in particular the linear veloc- 
ity variance can show IR divergences for n < — 1, which 
cancel out for equal-time statistics (Vishniac 1983; Jain & 
Bertschinger 1996). This suggests that it would make more 
sense to factor out a Gaussian factor of the form 



ti{k;D^,D2) 



47rfc3 



n=l 



H„(fc;7?i,Z?2) 



,(187) 



which is equal to unity at equal times. This also agrees with 
the functional form of the exact two-point function (|112p 
and of the simple approximation (|199p below. We show the 
results of the expansion p87p in Fig. [TO] again using the ex- 
act response function (jl33p . The high-order terms now grow 
as power laws at high k for the equal-time power (clearly 
they would still decay as a Gaussian for unequal times) , but 
the expansion does not fare better than the straightforward 
expansion obtained from Eq. (|179p displayed in Fig. [TOl 
This shows that using a reasonable Gaussian decay ansatz 
(which must disappear at equal times) is not sufficient to 
bring a significant improvement over previous expansions. 



9. High /c limit ? 

To improve the behavior of expansion schemes, Crocce & 
Scoccimarro (2006b) suggest using a response function that 
matches both the low-i behavior (obtained by perturbative 
expansions) and the Gaussian decay at high t. Indeed, for 
the case of the gravitational dynamics where the exact re- 
sponse R is not known, one can still obtain the low-t be- 
havior by perturbative expansions, such as those described 
above, by summing over higher-loop diagrams; see Crocce 
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& Scoccimarro (2006a) and Valageas (2007). On the other 
hand, Crocce & Scoccimarro (2006b) manage to resum a 
subset of diagrams in the high-fc limit, which gives rise to 
the same Gaussian decay e~* 1"^ as for the Zeldovich dy- 
namics. Then, one may hope that by using a "fit" for the 
response function that closely follows the expected behav- 
ior, one has done half the work needed to obtain a good 
prescription for the matter power spectrum, and all that is 
left is to use a good recipe for the self-energy H. In fact, we 
have shown in Sect. 18.31 and Fig. [TO] that this is not so sim- 
ple, because even using the exact response i? is not enough 
to improve the predictions for the power spectrum if we use 
a simple expansion over powers of Plq for the self-energy 
n, as in Ea. (|179p . Moreover, as shown in Sect. 18.41 using 
such a Gaussian cutoff for 11 is not sufficient either. 




Fig. 13. The expansion of the nonlinear field V' over the 
hnear growing mode t^jl from Eq. (|5T|) . up to order The 
filled circles are the vertex Ks^ whereas the white circles 
are the linear input 'i/'L- The numbers are the multiplicity 
factor associated with each diagram. 



Nevertheless, this section revisits the approximate re- 
summation performed in Crocce & Scoccimarro (2006b) to 
clarify its meaning and its shortcomings. Let us start from 
the general integral equation of motion ([?T|) . This equation 
can be solved perturbatively as an expansion over powers 
of ^/li of the form 



3„;,4 



(188) 



which can also be written in a diagrammatic manner, see 
Crocce & Scoccimarro (2006a) and Valageas (2001), as 
shown in Fig. 1131 In the high-fc limit, one assumes that 
all wavenumbers Wi associated with the linear fields '0l are 
much smaller than k except for one field (because of the 
conservation of momentum associated with the Dirac fac- 
tor (5D(ki -I- k2 - k) in the vertex K^, see Eq.j^Dl))- This 
assumes that the power is generated over some finite range 
of wavenumbers (i.e. the dynamics is not governed by an 
extended UV tail). Crocce & Scoccimarro (2006b) assume 
that the dominant contribution is provided by the diagrams 
shown in FigfTil where all low-Wi fields ^Li'^i) are directly 
connected to the "principal path" that joins the only one 
high-fc field to the root of the diagram. The idea is 
that in such diagrams the "principal path" only interacts 
with other fields ipLiyfi), which are pure linear growing 
modes, as they have not already suffered nonlinear interac- 
tions through the coupling vertex Ks- Then, these diagrams 



should maximize the cross-correlations since nonlinear in- 
teractions are expected to erase the memory of initial con- 
ditions. Therefore, the response function R may be domi- 
nated by these diagrams in the high-fc limit. 



O 







k k' 











Fig. 14. The diagrams assumed to dominate in the high-fc 
limit. We have fc' ~ fc, and all intermediate wavenumbers 
Wi are much smaller than fc. 



On the other hand, the diagrams of Fig[T4]are actually 
generated by the equation of motion 



2Ks{x; xi, X2).iJL {xi)il){x2) 



(189) 



where 'ijjf^{xi) is the linear growing mode restricted to low 
wavenumbers wi <C fc (the factor 2 comes from the fact that 
we can associate "ipL to either xi or X2 in Eq. (j51[) ). and we 
noted by a hat the approximate field ^ obtained with this 
high-fc limit. Note that and 'ijj^{xi) are independent 

Gaussian fields (since wi ^ fc) and Eq. (jl89p is now a linear 
equation for the field -0. Indeed, we can check that, by solv- 
ing Eq. (|189p as a perturbative series over powers of V'l, we 
recover the diagrams of Fig |14l In other words, by keeping 
only the diagrams of Fig[T3]we have actually approximated 
the nonlinear equation of motion (|5ip by the linear equa- 
tion of motion (|189p . Next, following Crocce & Scoccimarro 
(2006b) we note that the vertices 7"* of Eqs.(I2T])-(1211) satisfy 
at the high-fc limit, 



-» 1 : I]7';,„j(w, k)V'<„(w, -n' 



k.w 

'2^ 



e''(5Lo(w),(190) 



where we have only kept the terms of order k/w and ne- 
glected terms of order l,w/k, ... At this level, we can also 
replace the Dirac factor (5zj(w -I- k' — k) by SoCk' — k) in 
the vertex Kg, and using the second Eq. ([5^ and Eq. d^O]) . 
we can write Eq. (jl89p as 

/■ kw r - 

Sik, 77) = e''<5Lo(k)+ / dw— ,5Lo(w)e'' / d77',5(k, 77'), (191) 

J ^ J Tjl 

where we use the fact that ipL is also the linear growing 
mode. Here we consider the case where the initial conditions 
are set up at a finite time "qi as in Sect. 13.2.21 This linear 
equation can be solved through the expansion 

00 p 



— ' ft. 1 

5(k, ^) = €"5^0(1^) 4- e^Sm{k) YA\ \ -^'5lo(w,) 

(192) 



f) nix 

X / dr/ie''^ / d?72e''^. 



dr/pe''^ 



which can be resummed as: 

<5(k, D)=D (5Lo(k) pj^-^^) / dw^5^o(w) ^ 



(193) 



where we use the time-coordinate D — e^. Let us recall 
that SLo(k) and ^lo(w) must be treated as independent 
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Gaussian variables in Ea. (|193p . Then, if we define a re- 
sponse i?/(k, rj; k') in a fashion similar to Eg. ([55]) by 



i?/(k,r;;k')- 



where V is the functional derivative, we obtain 



(194) 



R 



J = 5^(k-k') — (e^^-^^)/''^^'-^^)) 
Di 

5^(k-k')— e-^(^-^^)''=''^'. 
^ ' Di 



(195) 



Thus we recover the exact Gaussian decay at high k ob- 
tained by Crocce & Scoccimarro (2006b). Note that since 
we have restricted the initial conditions to the linear grow- 
ing mode in Eq. (|19ip , the response Rj of Eq. p94p actually 
corresponds to the sum TZu +TZi2 of the components of the 
response defined in Eq. ([63l) . The advantage of this formu- 
lation is that one may see more clearly through Eqs. p89p - 
P93p the meaning of the assumptions involved in this high- 
k limit. In particular, it is interesting to compare Eq. (|193p 
with the exact result (|98p . which reads as 



S(k,D) 



e-»k.q D J dv 



■ i5i,o(w) 



i2n) 



(196) 



Then, if we assume that we can split the integral over w 
into low and high wavenumber parts w < A and w > A, 
such that most of the power is associated with w < A and 
the high-wavenumber contribution is small, we can expand 
the exponential over this high-wavenumber part as 



S{k,D) 



p-.k.q D J dw e™- il^ 5^„(w) 

(27r)3 



r 1 ' 

l + D dw'e™'<i-^Jio(w') 

Jiu'>A W 



(197) 



Next, if we assume that this expression is dominated by 
q ^ 1/k, we can neglect the factor e^^'^ in the exponent in 
the high-fc limit fc 3> A. Then the integration over q yields 
the Dirac factor S o [W — k) to obtain at lowest order (the 
factor 1 does not contribute) 



(198) 



Thus we recover Eq. (|193p with Dj and letting A — > oo . 
Of course, the assumptions involved in the derivation of 
Eq. (|198p from Ea. (|196p are identical to those involved in 
the derivation of Eq. (|193p . To check whether these assump- 
tions are valid, we can compare the nonlinear power spec- 
trum predicted by Eq. (jl93p with the exact power studied 
in Sect.O This gives (with Dj = 0) 



A'{k;Di,D2) = Ai e-5 



(199) 



We again recover the exact Gaussian decay at high k for un- 
equal times, which corresponds to the exponential term in 
the exact expression (|112p : but for equal times, we merely 
get back the linear power A|^. In fact, the analysis of 
Sect. 16.21 shows that the assumptions underlying Eq. (|198p 
are not valid. For a power-law linear power spectrum (|114p . 
the derivation leading to Eq. (|117p shows that, in the highly 
nonlinear regime the nonlinear power at wavenumber k is 
associated with scales q ~ i^-i+{n+3)/(n+i) ^^-^^ with linear 



wavenumbers w ^ fci-(n+3)/(n+i)^ Thus, for -3 < n < -1, 
where the system is well-defined, we find that the power at 
a nonlinear wavenumber k is generated by linear wavenum- 
bers w, which actually grow fastciQ than k instead of being 
restricted to a finite range w < A. Therefore, one cannot 
define a fixed cutoff A beyond which nonlinear interactions 
are negligible so that we can expand the high-w part as in 
Eq. p97p . The latter integral is actually large in the nonlin- 
ear regime, and one needs to take its full non-perturbative 
expression into account to obtain the non-trivial scaling of 
Eq.miHl). 

As was already clear from the linear Eqs. (|189p . (|19ip . 
the high-fc approximations described above neglect the 
nonlinear interactions associated with high wavenumbers, 
which actually govern the formation of large-scale struc- 
tures on small scales, and only keep track of the overall 
displacement associated with the large-scale velocity field. 
This is why the low-fc nonlinear interactions could be re- 
summed as e^^-'^i^^^'*'' °'" , which only involves the variance 
ay of the linear displacement field. This is clearly an impor- 
tant feature of the dynamics for unequal times, especially 
for the simple Zeldovich dynamics studied in this paper 
where the exact trajectories satisfy the simple law ([TU)) and 
the response function p33p happens to be fully determined 
by a similar advection process. In fact, Eq. p93p reads in 
real space (with Dj = 0) as 



,5(x,i?) = 5L(x-Si(q = 0,i?)), 



(200) 



where s^, = D+s^o is the linear displacement field of 
Eq. pH)) . As expected, Eq. (|200p explicitly shows that the 
approximation (5d(w + k' — k) ~ (5£)(k' — k) in the vertex 
Kg used to obtain Eq. (|19ip and to simplify the diagrams 
of Fig. [T3] has broken the invariance through translations 
of the system. As seen above, this invariance is restored as 
in Eqs. (|195p and (|199p by treating Si(q = 0) as an inde- 
pendent Gaussian random variable. Equation (j200p clearly 
shows that the effective dynamics associated with these 
high-fc approximations is simply the uniform advection of 
the linear density field by the linear velocity at q = 0. After 
averaging over the Gaussian initial conditions, this random 
displacement of large-scale structures leads to an apparent 

,2 2 

"diffusion" in the form of a Gaussian decay e~ for both 
the response function and the different-time correlation, see 
Eqs. (|195p and (|199p . This apparent loss of memory happens 
to dominate the different-time behavior of two-point func- 
tions for the Zeldovich dynamics, as seen from Eqs. (|112p 



^ That the power at a nonlinear wavenumber k mostly comes 
from higher linear wavenumbers as shown by the explicit 



expression (|109|l analyzed in Sect. 16.21 may seem a bit counter- 
intuitive as one may expect a "direct cascade" from larger to 
smaller scales. However, it appears that the actual process is 
somewhat more complicated within the highly nonlinear regime. 
In particular, as analyzed in Taylor & Hamilton (1996) (see also 
Schneider & Bartelmann 1995), the equal-time nonlinear power 
behaves as A^(fc,D) oc D~'^ (which is independent of fc) in the 
highly nonlinear regime if the spectral index is less than —3 on 
small scales (e.g. the linear power shows a high-fc cutoff Phik) oc 
fc^e^*" 1'^"). This scaling can be read from Ea. p05|l with m fixed 
and g ~ 1/ (-Dfc) , which shows that the power comes from a fixed 
range of wavenumbers w (e.g. w < K) associated with caustics 
(Schneider & Bartelmann 1995). However, for a spectral index 
larger than —3 on small scales the scaling is quite different as 
shown by Ea. (|118|l (see also Taylor & Hamilton 1996), and one 
cannot neglect smaller scales to describe the nonlinear evolution. 
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and (|133p . but it is actually disconnected from the building 
of matter clustering as is obvious from Eq. (|200p . 

It is clear that all steps going from Eqs. (|188p to p95p . 
and Eas. (|199p - (|200|) can be applied identically to the exact 
gravitational dynamics, see Crocce & Scoccimarro (2006b). 
As seen above, the high-fc approximations associated with 
Eq. (jl93p cannot aim at capturing the physics of gravita- 
tional clustering but only the "diffusion" associated with 
the linear velocity variance, which affects different-time 
statistics. However, it is not obvious that the approximate 
response (I195P should again agree with the exact response 
in the high-fc limit, which may thus depart from such a 
Gaussian decay. The small-scale gravitational dynamics is 
indeed quite different from the simple Zeldovich dynam- 
ics and other processes may come into play. Nevertheless, 
the apparent loss of memory due to the almost uniform ran- 
dom displacement of large-scale structures by low-fc modes, 
which is captured by the simple dynamics (|200p . clearly ap- 
plies to the exact gravitational dynamics as well. Therefore, 
one can expect again a decay as fast as in Eq. (|195p al- 
though only wavelengths which are still linear would con- 
tribute to (7„ (it is not clear whether nonlinear wavelengths 
would give rise to a stronger or weaker decay as compared 
with Ea. (|195p ). On the other hand, for the exact gravi- 
tational dynamics the fluid equations break down beyond 
shell-crossing so that the small-scale limit associated with 
these hydrodynamical equations is not so well defined. 



10. 2PI effective action metliod 

We now investigate the 2PI effective action method pre- 
sented in Sect. 14.21 and described in detail in Valageas 
(2007). Thus, we need to solve the system of coupled equa- 
tions (17ni)-(1721) and dUD-dHH). Thanks to causahty, which 
leads to the Heaviside factor ^(771 — 772) within both R and 
S, we solve this system by moving forward over time. We 
refer the reader to Valageas (2007) for a description of the 
numerical scheme. We consider the ACDM cosmology asso- 
ciated with Eq. (|153p and we only investigate the one- loop 
predictions. 

We first display in Fig. [12] the evolution forward over 
time Di of the response R{k; Di, 1)2). We can see that the 
nonlinear response exhibits oscillations as for the steepest- 
descent result (|163l) but its amplitude now decays as an in- 
verse power-law at large times Di instead of following the 
linear envelope (at one-loop order). As shown in Sect. 6.1 
of Valageas (2007) this behavior is due to the nonlinear- 
ity of the Schwinger-Dyson equation for the response R. 
Of course, in the weakly nonlinear regime we also recover 
the exact response (|133p (dashed lines). Next, we display in 
Fig. [in] the response function as a function of wavenumber 
k. In agreement with Fig. llSi we again obtain damped oscil- 
lations in the nonlinear regime. This is a clear improvement 
over both the standard perturbative expansion displayed in 
Fig. [3] and the steepest-descent result displayed in Fig. [6] 
This behavior is identical to the one obtained for the grav- 
itational dynamics studied in Valageas (2007). 

Finally, we show the logarithmic power A^{k; Di, D2) 
as a function of wavenumber k in Figs. [iTl and [TSl We com- 
pare the 2PI effective action prediction at one-loop order 
with the linear power, the exact nonlinear power obtained 
from Eq. (|109p , and the usual one-loop result obtained from 
standard perturbative analysis. The latter may also be ob- 
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Fig. 15. The nonlinear response R{k; Di, D2) (solid lines) 
as a function of forward time Di, for D2 — 0.32 (i.e. Z2 — 3) 
and wavenumbers k = 0.3 (left panel) and 3 x ft Mpc~^ 
(right panel) for the 2PI effective action method at one- 
loop order. For comparison we also plot the exact response 
i?NL (dashed lines). 
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Fig. 16. The nonlinear response function i?(fc; Di, D2) 
(solid lines) as a function of wavenumber fc, at times zi = 
0, Z2 = 3, for the 2PI effective action method at one-loop or- 
der. We also plot the exact nonlinear response i?NL (dashed 
lines). 

tained by expanding Eq. (|109p up to order P|q, and it reads 
for the Zeldovich dynamics as, 

P^'°°nk;D,,D2)=PL+P22 + Pi3, (201) 
with 

Pl = DiD2PLoik), (202) 
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on small scales and for large time separations; however, it 
is only a power-law decay instead of the exact Gaussian 
damping seen in Eq. (|112p . On the other hand, for equal 
times Di — D2 = D, we obtain a steady growth of the 
power A^(fc) in between the linear prediction A| and the 
usual one-loop prediction A^j^^p, see Fig.[T7l This is due to 
the contributions of nearby times D'l ~ D'2 ~ D in the last 
term of Eq. (|83p . which are not damped because of their 
small time-difference. Note that the cancellation at equal 
times of the damping associated with the decay of the re- 
sponse function is qualitatively correct, as shown from the 
exact nonlinear solution studied in sects. \6A\ and 16. 2[ see 
Eq. (|112p . Therefore, at one- loop order, the 2PI effective 
action method shows a significant qualitative improvement 
over the standard perturbative expansions of Sect. [7] and 
the direct steepest-descent method. However, it does not 
manage to predict the high-fc smooth power-law decay of 
the equal-time power A^(fc). 



Fig. 17. The logarithmic power A^(fc) (solid line) at red- 
shift z = 0, that is, at equal times zi = Z2 = 0. We also 
display the linear power A| (dotted line), the usual one- 
loop perturbative result A^j^^p of Eq. (|20ip (dotted hue) 
and the exact nonlinear power A^^ oi Ea. p09p (dashed 
line) . 
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Fig. 18. The logarithmic power A^(fc) at unequal times 
{Di = 1, D2 = 0.32) (i.e. zi = 0, Z2 = 3). 

Pi3 = - ^'^^' D,D2PLo{k)k'al (203) 

The results match the steepest-descent predictions, as well 
as the usual one- loop power (|20ip . on large scales. On small 
scales, contrary to the usual one-loop power (|20ip and the 
steepest-descent predictions, the 2PI effective action meth- 
ods yields a logarithmic power A^(fc; Di, D2), which decays 
for different times Di D2, see Fig.[THl This high-fc power- 
law damping is due to the decay of the response function 
already shown in Figs. ll5l and ll6l leading to a decorrelation 



11. Simple nonlinear schemes associated with the 
2PI effective action method 

11.1. Response function 
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Re(r) Re(r) 

Fig. 19. Left panel: the trajectories of the four roots of 
Eq. (|212p as we follow s along the real axis from s = +00 
down to s = 0. Right panel: the trajectories of the six roots 
of Ea. (|213p as we follow s along the line Im(s) = 1.052 from 
Re(s) = -|-oo down to Re(s) = 0. The root of interest start- 
ing from f = "collides" with a second root and changes 
direction by 7r/2. This corresponds to a singularity for the 
implicit fmiction f{s). 

For the 2PI effective action method, it is not straightfor- 
ward to obtain the self-energy E at a given order from the 
exact expressions (I140p or (|142p . Indeed, for the steepest- 
descent approach investigated in Sect. [51 the self-energy at 
a given order simply corresponds to the truncation of its ex- 
pansion over powers of Plo', therefore, it could be directly 
obtained by expanding the exact result. By contrast, within 
the 2PI effective action scheme, the self-energy is obtained 
from a diagrammatic expansion in terms of the nonlinear 
response R and correlation G (defined self-consistently at 
this order). Then, the exact expressions of the response 
R and correlation G are not sufficient to fully define the 
equations associated with the 2PI effective action method 
at any order, so one must go back to its diagrammatic def- 
inition. To avoid this complication, and to take advantage 
of the known expressions of the exact two-point functions, 
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which allowed us to bypass the computation of high-order 
diagrams in Sect.[51 we investigate here a nonlinear expan- 
sion that is not identical to the 2P1 effective action but 
is expected to show a similar behavior. Thus, we look for 
an expansion of the self-energy S over the nonlinear re- 
sponse R. A simple way to build such an expansion is to 
use the expansions over 1/s of the Laplace transforms f(s) 
and (t(s). Thus, from Eq. (|164p and the Laplace transform 
of Eg. (fits)) , we have 
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(205) 
(206) 



Then, the series (j205l) may be inverted as 

- = f + f^ +3f'' ~20f^ + ... (207) 
s 



Composing this expansion with Eq. (|206p . we obtain 

a = f -r^ +Af^ -27f^ + 2'i8r^ + ... (208) 

This provides an expansion of the self-energy S in terms 
of the nonlinear response R. In real space this yields 
multiple integrals over r(t), 

a{t)=r{t)-f dh f dt2r{t-ti)r{ti-t2)r{t2) + ..,{209) 
Jo Jo 

in a fashion similar to what would be obtained for the dia- 
grammatic expansion associated with the 2PI effective ac- 
tion. The main difference is that the expansion (|208p only 
involves the response R, whereas the 2PI effective action 
expansion involves both R and G, as in Ea. ([5T|) . Note that 
this shows that one can define several nonlinear expansion 
schemes. That it is possible to write a simple expansion 
such as (|208p is due to the exact response R only depending 
on the linear power spectrum through the velocity disper- 
sion a^. This is not the case for the gravitational dynamics 
where it may not be possible to write an expansion for S 
only in terms of R. Next, truncating the expansion (|208p at 
a given order and substituting it into Ea. (jl36p . we obtain 
at order p = 1: 



r, 



1 = 0, 



+ 4 - 



However, we now find that the solutions f'^^ (s) de- 
fined from these equations have singularities in the right- 
hand half-plane Re(s) > 0. For order p = 2, this may 
be directly seen from the explicit solution of Eq. (|212p or 
from the behavior of the four roots ri{s) as a function 
of s shown in left panel of Fig. [191 Indeed, as we follow 
s along the real axis from +oo to 0, the root of interest 
fi(s) that starts from fi = at s = -|-oo "collides" with 
a second root r2 at ~ 0.87 (s* ~ 0.94) and afterwards 
forms a pair of complex conjugates with f2- This is asso- 
ciated with a square-root singularity for r(s). (A simple 
example is provided by the polynomial — s = with 
a singularity at f* = 0, s* = 0.) For order p — '3, the 
right panel of Fig. where we follow s along the line 
Im(s) = 1.052 from Re(s) = -l-oo down to Re(s) = 0, shows 
that we again have a singularity at the complex points 
s* ~ 1.17 ± 1.05i, ~ 0.56 =F 0-31«. These singularities 
yield exponential factors e''** for the response r(t), which 
grow at large t. Therefore, although the expansion ()208p is 
much better than the steepest-descent approach (I163P at 
order p = I, since it exhibits a damping in the nonlinear 
regime, it shows as for expansion (|163p growing exponen- 
tials at higher orders. Thus, in this sense this nonlinear 
expansion is not well-behaved. We can expect that a simi- 
lar problem occurs for the 2PI effective action approach at 
high orders. 

11.2. Correlation function 

We now investigate nonlinear schemes for the two-point 
correlation G. As for the response function, we look for a 
simple nonlinear expansion that bypasses the need to com- 
pute high-order diagrams but that follows the structure of 
the 2PI effective action method. This is not as straightfor- 
ward as for the response R because the two-point corre- 
lation G cannot be written in terms of a one-dimensional 
function such as r{t) for the response R. Thus, we focus on 
the equal-time nonlinear power for the case of a power-law 
linear power spectrum and we write 



A^{k;D) ^ Alg{t) with t^DjA^ 



LO 



(214) 



(210) 



which defines the time-variable t used in this section and 
the function g{t). From Eas. (|156p and (|175p . we have for 

n = -2 



The root of the polynomial of degree 2p which must be 
chosen, is the one that is consistent with the expansion 
(|205p for s — > oo. Equation (|210p is a well-known Laplace 
transform, and wc obtain 



,(1) 



it) 



Ji{2t) 



t 



(211) 



Note that this expression is also obtained as a simple ap- 
proximation for the one-loop 2PI effective action approach 
for the gravitational dynamics in Valageas (2007). We also 
recover the damped oscillations obtained within the 2PI ef- 
fective action method displayed in Figs. [T51 and [TBI Thus, 
as expected this nonlinear expansion and the 2PI effective 
action expansion show the same behavior at order p — \. 
At orders p = 2, 3, we obtain the polynomial equations: 



p^2 
p = 3 



f2 
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sf - 1 = 
+ sf - 1 = 0. 



(212) 
(213) 
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which yields for the Laplace transform defined as in 
Eq.([n7D: 



5(s)=^^(2p-2)!.-2^+i = i 

, p! s 

p=i 

This series can be inverted as: 



3^ 
s 32s3 
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s 32 ^ 1024 

Next, the derivative of g{t) verifies 

oo 

.9'w-E^(2p-2)^2^-^ 

p=2 P- 



4s5 



3tt^ o 3(2567r2 + 9tt*) . 
9 + TT^. -9^ 



..(216) 
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(218) 
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and the Laplace transform of this equation reads 



sg{s) - 1 
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pi 

p=2 f 

32s2 " 



f (2p-2)!s-2p+2 
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512s6 



(219) 



which could also be obtained from Eq. (|216p . Then, substi- 
tuting the series (PT7|) into Eq. (PTg|) gives 



3^2 



37r2(128 + 37r2)^4 
512 ^ 



(220) 



Thus, as for the response function studied in Sect. 111.11 
we have obtained a simple nonlinear expansion scheme for 
the two-point correlation G (i.e. for the nonlinear power 
A^(fc; D)). However, since the functional dependence of the 
two-point correlation G is not as simple as for response R, 
the expansion (|220p only applies to the equal-time power, 
and it is not built from the self-energy H. Nevertheless, 
going back to real space, Eq. (|220p yields an integro- 
differential equation for G such as Eg. ([70]) . with a fixed 
differential term on the left hand side and a series of mul- 
tiple integrals over g{t) on the right hand side. At order 
p = 0, the right hand side is zero and we recover the linear 
solution g^^'{t) = 1. At order p = 1, we obtain 



which gives 



37r2 



3TTt 



2 2 



(221) 



(222) 



where Ji is the modified Bessel function of order 1. Thus, 
the nonlinear power A^^^^ shows an exponential growth in 
the highly nonlinear regime. We can check that at order 
p = 2 we again have an exponential growth (with oscil- 
lations since Im(s*) ^ where s* is the location of the 
singularity of ^'■^''(s)). Therefore, the nonlinear expansion 
l|220p is not better behaved than the "linear" expansions 
associated with the steepest-descent approach studied in 
Sect. 



12. Simple nonlinear schemes associated with the 
running with the high-Zc cutoff 

12.1. Response function 

In a fashion similar to Sect. 1111 we now investigate some 
nonlinear schemes that may be built in the spirit of the 
method outlined in Sect. O where we considered the de- 
pendence of the system on a high-fc cutoff A. In order to 
separate the dependence on A, we now write the response 
function as 



R^RLr{t,uj^) with t ^ Di - D2 



■J Jo 



(223) 



(224) 



Thus, the exact response function r{t, lo^) is, from Eq. (|133p . 
r(i,t^2) ^ g-c^^tV2_ (225) 



Defining the Laplace transform with respect to t as in 
Eq.ldnZD, we obtain 



f (s, ^2) = ^(-1)P {2p - 1)!! w^p 5-2p-i^ 
while the derivative with respect to A is 



(226) 



{-If (2p- l)!!pw2(p-i), 



-2p-l 



(227) 



Inverting the series (|226p and substituting into Eq. (l227p 
gives the expansion 



df 



d^ 
dA 
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(228) 



At order p = 0, we have df/dA = and we recover the 
linear response (since we impose r(A = 0) — r^). At order 
p = 1 we obtain 



df 
dX 



d^ 
dA 



(229) 



which is similar to Eq. ((9T|) once we go back to real t- 
space, which leads to a double integral over time (as in 
Eq. pIS)) ). We can recover Eq.^ by noting that the hn- 
ear response is ^^(s) = 1/s, whence f| = l/2d^fL/ds^. 
Substituting this result into Eq. (|229p gives dr/dA = 
— l/2(da;^/dA)d^fi/ds^. The inverse Laplace transform of 
this equation gives back Eq. f^^ . It is clear that this pro- 
cedure is not systematic and only applies to the linear re- 
sponse. This is why we could reduce the three response 
functions on the right hand side of Eq.(PT|) to the one re- 
sponse function on the right hand side of Eq. ([g^ . Then, 
replacing Rl by i? in Eq. l^^ is clearly correct at the order 
at which the right hand side of Eq. ipO]) was truncated, 
but this method cannot be extended to higher orders in 
a systematic fashion. Thus, let us consider Eq. (|229p with 
keeping the right hand side as it comes out from the system- 
atic expansion obtained in (|228p . Integrating over A gives 



1 



Vs2 + 2u;2 



whence r'^^\t) = Jo(\/2cjt). (230) 



Thus, we obtain decaying oscillations into the nonlinear 
regime, which are damped as [k{Di — _D2)]~^^^- Note that 
the damping is smaller than for the 2PI effective action 
method where the simplified expansion (|208p gave Eq. (|21ip , 
which decays as [k{Di — D2)]~^^'^. At next order p = 2 we 
obtain 

df duj^, n ^ , clf 



-f^ 



3a;2p. (231) 



Then, looking for a solution of the form 

r{t,u;^) = p{Lut), f{s,iu^) = -p{y) with y = -, (232) 

we obtain for the Laplace transform p{y) the equation 
p + yp' = 2p^ -6p^. (233) 



Note that Eq. (|233p no longer involves an explicit depen- 
dence on to. It can be solved in implicit form as 

y ^ i(i_2/5V6/54)i/4eH*^"^~"'"'"(T^)]^^^^^(234) 
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Going back to s 
the point 



ujy, we see that f(s,w^) is singular at 13. Weakly nonlinear scales 



with 



y* 



where |f| = oo. Therefore, the response function obtained 
at order p — 2 grows into the nonUnear regime as r*-^) ^ 
6"'^ which gives i?^^) ^ gjy.fc ^„(Di-C2)^ Thus, as for the 
nonhnear scheme of Sect . 1 1 1 . 1] associated with the 2PI effec- 
tive action method, we find that at high orders, the expan- 
sion (|228p gives rise to response functions that exhibit an 
exponential growth into the nonlinear regime, even though 
at lowest-order p = 1, it managed to provide a nonlinear 
damping. Therefore, that the Gaussian decay could be re- 
covered from Eq. (|5^ is not due to the good convergence 
properties of the method outlined in Sect. [5] As discussed 
in Sect. [5l from Eq. (|94|) it merely comes from the succes- 
sive steps that have been performed using the properties 
of the linear response until one gets the linear Eg . ([92| . 
which is correct at lowest order and which has the desired 
solution. However, these intermediate steps cannot be di- 
rectly extended to high orders without ambiguities, and the 
systematic nonlinear expansion (|228p does not recover this 
damping at high orders. 

12.2. Correlation function 

For the correlation function G, following the procedure of 
Sects. [TL^ and [Tm we write 



A'^{k;D) = Al g{t, Ala) with t = D. 



(236) 



Then, again introducing the Laplace transform with respect 
to t, inverting the series g{s, A|^q) — l/s + .. and substitut- 
ing into dg/dA gives the expansion: 



dA 



dA|o 
dA 



3^ 
32 



~9' 



37r2(512 + 97r2) 
1024 



(237) 



At order p = we have the linear correlation dg/dA — 0, 
g — 1/s and g{t) = 1. At order p = 1 we obtain 



In the previous sections, we have investigated the conver- 
gence properties of several expansion schemes that may 
be used to study gravitational clustering in the expand- 
ing Universe, applied to the case of the Zeldovich dynam- 
ics where exact results can be obtained. In this section, 
we complete this study by a brief description of the re- 
sults obtained at one-loop order on weakly nonlinear scales. 
Indeed, an accurate prediction for the matter power spec- 
trum on these scales is of great practical interest for sev- 
eral cosmological probes, such as baryonic acoustic oscilla- 
tions (Eisenstein et al. 1998, 2005) and weak-lensing shear 
(Munshi et al. 2007). For the exact gravitational dynam- 
ics, the various expansion schemes must be compared with 
numerical simulations, which is not very convenient to eval- 
uate their power. Therefore, it is interesting to check the ac- 
curacy and the behavior of these expansion methods against 
the Zeldovich dynamics. 

13.1. Linear expansion schemes 




dg 



32 

which gives 



3^ ~9^'H-^) 



(238) 



ffW(t,Aio)-/o 



(239) 



where Iq is the modified Bessel function of order 0. Thus, 
as for the simple nonlinear scheme associated with the 
2PI effective action method of Sect. I11.2[ we obtain at or- 
der p = 1 a nonlinear equal-time power A^ which shows 
an exponential growth in the nonlinear regime. Following 
the method leading to Ea. (|234p . we can actually integrate 
the nonlinear equation (j237p at any order, using the scal- 
ing g{s, A|q) = 7(s/Alo)/Alo, which transforms Ea. (|237p 
into an implicit equation giving y = s/A^q as the integral 
of a rational function of 7. We can check that, at order 
p = 2, we again have an exponential growth into the non- 
linear regime (with oscillations because lm{y^) ^ 0). 



Fig. 20. The power spectrum P{k) divided by a smooth 
linear power p™°°*h a,t redshift z = 0. We display the lin- 
ear power Pbik) (dashed line), the exact nonlinear power 
Pnl of Eq. (|109p . the standard 1-loop result (dotted line) 
of Eq. (|20ip , and the steepest-descent result of Eq. ((83l) (up- 
per solid line "s.d."). We also show the results obtained 
by adding a Gaussian factor to the standard perturbative 
result as in Ea. (|159[) (lower dot-dashed line) or by using 
the exact nonlinear response -Rnl in Ea. (j83p within the 
steepest-descent scheme (upper dashed line), as in Fig. [TUl 



We first consider "linear" expansion schemes, that is, 
methods that give rise to expansions in terms of the linear 
power spectrum, such as the standard perturbative expan- 
sions of Sect. Hand the direct steepest-descent methods of 
Sect.m We focus on the equal-time power spectrum for the 
ACDM universe described in the first paragraph of Sect. [71 
We used the CAMB Boltzmann code (Lewis et al. 2000) to 
obtain the linear power spectrum with the baryonic acous- 
tic oscillations. In order to magnify the difference between 
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Fig. 21. The power spectrum P{k) divided by a smooth 
hnear power p™o°th as in Fig.[2Ql but at redshift z = 1. 

various schemes, we show in Figs. and HT] the nonhn- 
ear power divided by the linear power p£™ooth associated 
to a smooth power spectrum without baryonic oscillations, 
taken from Eisenstein & Hu (1998). We compared the re- 
sults of various expansion schemes with the exact nonlinear 
power (solid line labeled Pnl) obtained from the numeri- 
cal integration of Eq. (|109p for the ACDM power spectrum. 
First, we see in Fig. [20] that all schemes agree with the ex- 
act power up to /c ~ O.lh Mpc^^ at redshift z = and 
follow the departure from the linear power P^. On smaller 
scales, the various expansion schemes deviate from one an- 
other and from the exact nonlinear power. It is interest- 
ing to note that using the exact nonlinear response (|133p 
within the steepest-descent scheme does not improve the 
agreement over the original steepest-descent scheme where 
we use the response function predicted at the same order. 
On the other hand, it appears that the standard one-loop 
expansion provides the best results at this order. 

Following Ea. (|159|) and Fig.[5l we also consider the ex- 
pansion defined from the standard perturbative series by 
factorizing a Gaussian decay: 

PG:Z.ik) = e-^'-' [Pl + P,2 + Pi3 + DWPl] 

= e-^"^\PL+P22). (240) 

In agreement with Sect. [7] and Crocce & Scoccimarro 
(2006a), this gives a positive power whatever the shape of 
the linear power spectrum; however, Fig. l20l shows that this 
does not necessarily improve the accuracy as compared with 
the usual one- loop result (|20ip . We show the power spec- 
trum obtained at redshift z = 1 in Fig. [21] Of course the 
various expansion schemes agree more closely on the same 
scales with the exact result, since we are closer to the linear 
regime. We can see that we recover the same behaviors as 
in Fig. [20] 

The behavior of these various expansion schemes at 
higher orders was investigated in sects. [7][8] We found 
that no linear scheme provides a significant improvement 
over the standard perturbative expansion. Factorizing a 
Gaussian decaying term as in Eq. (|240p seemed to give 
a small improvement in Fig. [5] but as discussed below 



Fig. [5] this is not very robust and would not apply to any 
power spectrum. On the other hand, Fig. [TU] shows that 
even using the exact nonlinear response within the direct 
steepest-descent scheme does not generically provide a sig- 
nificant improvement either. Nevertheless, the most promis- 
ing scheme within this framework is probably to combine 
a good ansatz for the response function with such expan- 
sions (see also Sect. 113.^ below). Then, one could hope to 
gain from the higher orders, while the imposed decay of the 
response function could tame the increasingly fast growth 
at higher orders obtained in the standard perturbative ex- 
pansion. An application of such a strategy was presented 
in Crocce & Scoccimarro (2007) for a ACDM universe, and 
a significant improvement over the standard perturbative 
expansion was obtained up to two-loop order. 

13.2. Nonlinear expansion schemes 
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Fig. 22. The power spectrum P(fc) divided by a smooth 
linear power p™ooth redshift z = 0. We display the linear 
power PLik) (dashed line), the exact nonlinear power Pnl 
of Eq. fTUg)) . and the 2PI effective action result of Eq. (|55)) 
(upper solid line "2PI"). We also show the results obtained 
by using the exact nonlinear response i?NL in Eq. (jS^)) within 
the 2PI scheme (lower dashed line), or by using the linear 
two-point correlation to compute both S and 11 (upper dot- 
dashed line) or only E (dotted line). 

Finally, we study in this section "nonlinear" expansion 
schemes, that is, methods that give rise to expansions in 
terms of the nonlinear two-point functions R and G, such 
as the 2PI effective action method of Sects. 1421 and [TOl We 
show our results for the power spectrum (divided again by 
psmooth) redshifts z = 0, 1 in Figs. ^ and [H First, 
we note that the 2PI effective action result at one-loop or- 
der overestimates the power spectrum on weakly nonlin- 
ear scales. Both the standard one-loop result and the di- 
rect steepest-descent method actually work better in this 
range. This behavior can be traced back to the damping 
self-energy E. Indeed, since E cx RG at one- loop order (see 
Ea. (|8T|l ). the decay of the response and of the two-point 
correlation at high k (shown in Figs. [TBI and [T5)) ) leads to a 
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Fig. 23. The power spectrum P{k) divided by a smooth 
hnear power p£"io°th as in Fig. [221 but at redshift z = 1. 



smaller at high k as compared with the Sq obtained 
for the direct steepest-descent method. This in turns yields 
a response R that is somewhat larger than for the steepest- 
descent method in the weakly nonlinear regime (but smaller 
at high k where it decays) whence a two-point correlation 
G that is somewhat larger from Eg. ([83)) . Of course, this 
slight overestimate of R and underestimate of E will be cor- 
rected by higher-order terms. For the Zeldovich dynamics, 
we know that R and S actually depend only on , so that 
terms that depend on other integrals over PLoik) must can- 
cel out in the full resummation. However, this discrepancy 
makes the one-loop 2PI effective action result insufficient 
for practical purposes (at least at the one-loop order). 

We also display the results obtained when both the self- 
energies E and 11 are obtained from the linear correlation 
Go (so that only the response R coupled to E is obtained 
from nonlinear equations) and when only the self-energy 
E is obtained from Gq (so that the nonlinear systems for 
the pairs {E,i?} and {11, G} are decoupled). We see that 
these two methods give results very close to the original 
2PI effective action prediction where all two-point functions 
{E, n, R, G} are coupled. 

On the other hand, we also show the power spectrum 
obtained from the coupled equations (|82p and ([55)) , when we 
use the exact nonlinear response (|133p . We can see that this 
significantly improves the agreement with the exact power 
spectrum, and the result is slightly better than the steepest- 
descent prediction shown in Fig. [201 but it is still a bit less 
accurate than the standard perturbative result at z = 0. 
However, Fig. [221 shows that at z = 1 the power obtained 
in this fashion is slightly more accurate than the standard 
perturbative result. Therefore, contrary to the case of the 
steepest-descent method studied in Sect. ll3.Tl using the ex- 
act nonlinear response improves the prediction for the two- 
point correlation and provides a scheme that can be com- 
petitive with the usual perturbative expansion on weakly 
nonlinear scales. 

At higher orders, the nonlinear 2PI expansion scheme 
- especially if it uses the exact nonlinear response or a re- 
liable ansatz (as for the lower dashed lines of Figs. [^Dl[^5)) 



- may provide an even greater improvement as compared 
with the standard perturbative expansion. The analysis of 
Sect. rTT|[T2l showed that nonlinear schemes can give un- 
wanted exponential growths at higher orders for both the 
response function and the power spectrum. However, this 
does not rule out a good convergence on weakly nonlinear 
scales. Nevertheless, this behavior and Figs. [20l[23l suggest 
that the most promising scheme would be to combine a 
good ansatz for the response with such nonlinear expan- 
sions. As for the linear schemes of Sect. 113.11 one could 
hope to benefit from the higher orders, while the imposed 
decay of the response function could partly restrain their 
increasingly fast growth. This could make the improvement 
over the standard perturbative expansion even more signif- 
icant at higher orders, but such a study is left for future 
work. 



14. Conclusion 

In this article we have applied to the Zeldovich dynam- 
ics various expansion schemes that may also be used for 
the gravitational dynamics in the expanding Universe. We 
derived the path-integral formalism that describes this sys- 
tem, starting either from the diff'erential or the integral form 
of the equations of motion, and we obtained the relation- 
ship between the associated response functions. These re- 
sponse functions describe the response of the system to a 
small perturbation applied at any time and they also en- 
code the memory of initial conditions. Next, we briefly de- 
scribed how to build various expansion schemes from these 
path integrals, such as large- expansions or running with 
a high-fc cutoff. All these results apply almost identically to 
the case of the gravitational dynamics. 

Then, we have derived the exact nonlinear two-point 
functions associated with the Zeldovich dynamics, taking 
advantage of the well-known exact solution of the equa- 
tions of motion. Whereas the equal-time nonlinear power 
decays as a power law in the highly nonlinear regime, 
different-time two-point functions, such as the response 
R{k;ti,t2) or the correlation G{k;ti,t2), show a Gaussian 
decay e"*^ '^v{Di-d.2) ji ^ where is the variance of the lin- 
ear velocity. This damping is associated with the uniform 
random displacement of particles between different times by 
long wavelength modes. This leads to an effective decorre- 
lation but it is not directly related to the matter clustering. 
In particular, tr^ can be made very large, and even infinite, 
through low-fc divergences, without affecting the equal-time 
matter power spectrum. This also means that the matter 
power spectrum P(fc; t) may still be very close to linear on 
scale k even though two-point functions such as i?(fe; ti^ti) 
may have already shown large deviations from their linear 
values on the same scale at previous times < t\ < t. 
Therefore, departures from linearity of different two-point 
functions are not necessarily related. 

Next, we have studied the standard perturbative expan- 
sions for both the response R and the logarithmic power 
A^(fc; t) (X k'^P{k; t), focussing on a power-law linear power 
spectrum n = —2 for the latter. The usual perturbative 
expansion being equivalent to a Taylor expansion it cannot 
capture the different-time Gaussian decay of R and A^, 
nor the equal-time power-law decay of A^, since it gives 
polynomial approximations of increasing order. In the spirit 
of Crocce & Scoccimarro (2006a), we noticed that for the 
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power A^, reorganizing the perturbative series by factor- 
izing a simple Gaussian term ^ '^^ gives an expansion 
which looks better behaved, as all terms become positive 
and there is a Gaussian damping. However, this procedure 
only works for linear power-spectra such that the scales 
which govern and are close (as for typical ACDM 
power-spectra). 

Then, we have studied the steepest-descent method de- 
rived from a large- AT expansion. At first order p = 1, there 
is some improvement for the response function, which re- 
mains bounded (as a cosine) instead of growing as fc^. (The 
same behavior is obtained for the gravitational dynamics 
(Valageas 2007).) However, at higher orders the expansion 
worsens as it exhibits an exponential growth in the nonlin- 
ear regime. We showed that this could be cured by using 
Fade approximants, which remain bounded at all orders as 
a sum of cosines (they do not decay but after integration the 
oscillations should produce some effective damping). The 
power A^ displays an exponential growth at orders p > 2 for 
the plain steepest-descent method and a polynomial growth 
using the Fade approximants. As for the standard perturba- 
tive expansion, one can reorganize the series by factorizing 
a Gaussian term £"(''^1+^2'''' ^"/^ in the self-energy H that 
describes the generation of power by nonlinear interactions. 
This appears to give slightly better results than with the 
standard perturbative expansion, but this procedure obvi- 
ously suffers from the same restrictions. Alternatively, one 
can factorize a Gaussian term such as e^^^i^^a)*: '^1/2 into 
the self-energy H, to reproduce the fact that the Gaussian 
damping must disappear for equal-time statistics. Then, 
one again obtains a polynomial growth into the highly non- 
linear regime for the power A^(fc;t), because of the con- 
tribution of mode couplings at recent times ti ~ t2 — t. 
Therefore, none of these methods shows very satisfactory 
global convergence properties. 

Next, we have discussed a high-Zc resummation proposed 
by Crocce & Scoccimarro (2006b) to improve the behavior 
of such expansion schemes. We showed that the partial re- 
summation involved in this procedure is equivalent to ap- 
proximating the nonlinear equation of motion by a linear 
equation. Using the high-fc asymptotic of the coupling ker- 
nel, as in Crocce & Scoccimarro (2006b), we derived the 
explicit expression of the nonlinear density field S{x, t) as- 
sociated with these approximations. Then, both the non- 
linear response R and power A^ are equal to their linear 
counterpart multiplied by the same different-time Gaussian 
decay. Thus, these approximations manage to capture the 
different-time Gaussian decay (associated with the advec- 
tion by long wavelengths) but they fail to capture the equal- 
time properties of the system. A close analysis of this pro- 
cedure shows that the underlying assumptions arc not valid 
because one cannot define a high-fc limit in a simple man- 
ner. (For a linear power spectrum with — 3<n<— lat 
high fc the nonlinear power at wavenumber k is generated 
by the highly nonlinear couplings of modes k' ^ k instead 
of being produced by small wavenumbers restricted to some 
finite range fc' < A.) The same caveats should apply to the 
gravitational dynamics, although it is not totally obvious 
whether the Gaussian decay at different times is exact in 
this case (but this is not necessarily important for practical 
purposes). 

Then, we have turned to nonlinear schemes, that is, ex- 
pansions over powers of nonlinear two-point functions, such 



as the 2FI effective action method built from a large— A^ 
expansion. At one- loop order [p = 1), we obtain damped 
oscillations (with a power-law decay) for the two-time re- 
sponse function and power A^(fc; ti, 1,2), but the equal-time 
power still grows on small scales. (One obtains the same 
behavior for the gravitational dynamics (Valageas 2007)). 
Then, from the exact two-point functions we have built sim- 
ple nonlinear expansion schemes that are similar to the 2PI 
effective action expansion, but that can be more easily han- 
dled analytically. We recover damped oscillations for the re- 
sponse function at order p = 1 , but we find an exponential 
growth at higher orders; the equal-time power also shows an 
exponential growth. Next, we investigated simple nonlinear 
expansion schemes associated with the evolution of the sys- 
tem with a high-fc cutoff A. We find a similar behavior as 
we again obtain damped oscillations at order p = 1 and ex- 
ponential growth at higher orders for the response function; 
the equal-time power also exhibits an exponential growth. 
To recover the Gaussian decay from this expansion at order 
p = I, one must introduce some further approximations, as 
in Matarrese & Pietroni (2007a,b), which are correct at 
this order but which do not give a systematic procedure. 
Therefore, it appears that at high orders p > 2 nonlinear 
methods do not fare much better than linear schemes. 

Finally, we have studied the quantitative predictions of 
these various schemes at one-loop order on weakly nonlin- 
ear scales for the equal-time matter power spectrum. All 
linear expansions agree with the exact results on quasi- 
linear scales (fc < Q.lh Mpc^^ at z = 0) where there is 
already a deviation of 10% from linear theory. On smaller 
scales they depart from each other and the standard per- 
turbation theory actually works best for the case studied 
in this article. Moreover, factorizing a Gaussian damping 
factor or using the exact response function does not im- 
prove the predictions obtained within this framework. On 
the other hand, the nonlinear 2FI effective action method 
overestimates the nonlinear power spectrum on these scales 
because of the nonlinear feedback involved by the coupled 
system obeyed by the response R. However, using the ex- 
act nonlinear response within this method now improves 
the agreement with the exact result and can be compet- 
itive with the standard perturbation expansion (but this 
depends on the exact shape of the linear power spectrum 
as in our case it is slightly better at z = 1 but slightly worse 
at z = 0). 

Since the equations of motion associated with the grav- 
itational and the Zeldovich dynamics are very close we can 
expect these results to apply to the gravitational case. (This 
is the case at one-loop order as shown by the compari- 
son with Valageas (2007).) We have found that none of 
these schemes shows good global convergence properties at 
high orders. Indeed, they all lead to polynomial or expo- 
nential growth into the nonlinear regime for the response 
function, except for the use of Fade approximants which 
gives bounded response functions with fast oscillations in 
the highly nonlinear regime. On the other hand, nonlin- 
ear schemes manage to reproduce the damping at one-loop 
order p = I but fail at higher orders. Next, no scheme 
manages to recover the power-law damping of the nonlin- 
ear matter power spectrum. They either display increasing 
growth at higher order or a Gaussian decay which is also 
somewhat artificial. 

Nevertheless, an expansion may still be very useful on 
weakly nonlinear scales even if it converges badly (or even 
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diverges) on highly nonhnear scales. There, we found that 
the best methods seem to be the standard perturbation 
theory or a nonlinear expansion where one uses the exact 
nonlinear response (with the Gaussian decay). These re- 
sults are somewhat disappointing, since it appears to be 
difficult to build systematic expansion schemes that sig- 
nificantly improve over the standard expansion. One may 
still obtain some improvement as in Crocce & Scoccimarro 
(2007) or Matarrese & Pietroni (2007a,b), but this requires 
some additional ingredients, such as the use of an ansatz, 
that shows a Gaussian decay, for the response function and 
in some cases for other two-point functions as well. On the 
other hand, the use of several expansion schemes can be of 
interest by itself, since they should be accurate at least over 
the range where they all agree. This allows one to obtain 
an estimate of their range of validity without the need to 
perform numerical simulations. 

In order to make progress, it appears that it may be 
advantageous for observational purposes to be guided by 
the expected behavior of two-point functions and to com- 
bine such systematic expansions with reasonable ansatze 
(e.g. Crocce & Scoccimarro 2007). From a theoretical per- 
spective, one may also look for different approaches. For 
instance, one could try to work directly with the Vlasov 
equation (Valageas 2004). However, this would make the 
computations significantly more difficult, and it is not clear 
whether it is not more efficient to stick to the hydrodynam- 
ical approach and to simply compute higher order terms, 
especially if one is mostly interested in weakly nonlinear 
scales. On the other hand, one may consider simpler ef- 
fective dynamics that attempt to go beyond shell crossing 
(based for instance on a Schrocdingcr equation, Widrow & 
Kaiser 1993). 
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